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THIRD YEAR (FINAL) PROGRESS REPORT 

This report covers technical progress during the third year of the NASA Space Physics 
Theory contract “The Structure and Dynamics of the Solar Corona, NAS5-96081, between 
NASA and Science Applications International Corporation, and covers the period June 16, 1998 
to August 15, 1999. This is also the final report for this contract. Under this contract SAIC, the 
University of California, Irvine (UCI), and the Jet Propulsion Laboratory (JPL), have conducted 
research into theoretical modeling of active regions, the solar corona, and the inner heliosphere, 
using the MHD model. During the three-year duration of this contract we have published 49 
articles in the scientific literature. These publications are listed in Section 3 of this report. In 
the Appendix we have attached reprints of selected articles. 

In the following sections we summarize our progress during the third year of the contract. 
Full descriptions of our work can be found in the cited publications, a few of which are attached 
to this report. 

1. PROGRESS REPORT 

Predicting the Structure of the Solar Corona for the August 11, 1999 Eclipse 

We have used our polytropic 3D MHD model to predict the structure of the solar corona 
during the total solar eclipse that occurred on August 11, 1999. Our coronal prediction was 
posted prior to the eclipse on the World Wide Web: 

http : //haven . saic . com/corona /model ing .html 
This prediction is our fourth in a series (previous predictions were made for the total solar 
eclipses in 1995, 1997, and 1998, as detailed on our Web page). The distinguishing feature of 
this prediction is that it occurred at a time when we were approaching solar maximum, so that 
the complexity of the Sun is much greater than for our previous predictions. This prediction 
was presented at the conference “The Last Total Solar Eclipse of the Millennium,” which was 
held in Istanbul, Turkey, August 13-15, 1999. Z. Mikic also witnessed the eclipse in Elazig, 
Turkey. A manuscript detailing the comparison was submitted for publication in the 
proceedings of this conference, and will be included in a future progress report. The 
comparison of our prediction and actual eclipse images is reasonably good. Figure 1 shows 
images of the predicted coronal polarized brightness, as well as the magnetic field line structure, 
and a comparison of the prediction with an image of the corona taken by Fred Espenak in 
Turkey. This coronal simulation is also being used to study the magnetic field as we approach 
solar maximum, as well as to compare with SOHO, WIND, and Ulysses solar wind 
observations. A preprint of past comparisons of our MHD model with eclipse observations has 
been published (Mikic et al. 1999), and is included in the Appendix. 

Fast Magnetic Reconnection in a 2D “Rosette” Configuration 

In conjunction with Dr. Samuel Vainshtein of the University of Chicago, we are studying 
the 2D reconnection of magnetic fields in a “rosette” topology in search of fast magnetic 


Comparison of a 3D MHD Coronal Prediction with 
an Image of the 11 August 1999 Total Solar Eclipse 



Fred Espenak’s Composite Image (Turkey) 



Predicted Polarization Brightness (MHD Model) 



Predicted Magnetic Field Lines (MHD Model) 


Figure 1. Comparison between a composite eclipse image created from photographs taken by Fred 
Espenak in Lake Hazar, Turkey (top) with the predicted polarization brightness of the simulated 
solar corona from our 3D MHD model (middle). The projected magnetic field lines from the 
model are also shown (bottom). Terrestrial (geocentric) north is vertically upward. The eclipse 
image is copyrighted 1999 by Fred Espenak. 
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reconnection. This is a follow-on to the work reported in the second quarter of this year in 
which we investigated the formation of a tangential discontinuity in the ideal MHD model. Fast 
reconnection is defined as one whose rate is weakly dependent on the magnitude of the plasma 
resistivity. We have found that, if the plasma viscosity is kept small enough, fast reconnection 
may occur. A series of simulations with r/ = v= 2.5 x 10' 5 , 5 x 10' 5 , 10 x 10' 5 , and 20 x 10' 5 
(where rj and v are the plasma resistivity and viscosity, respectively) show that the reconnection 
rate is almost independent of the resistivity. We are continuing these promising simulations to 
confirm that the reconnection rate is indeed fast. If this. is confirmed, this work would have 
impact on the applicability of reconnection as an energy release mechanism for coronal heating 
and solar flare initiation. 

Figure 2 shows the reconnection of two current-carrying flux bundles for the case w'hen 
rj= v = 5 x 10' 5 . Note that the two flux bundles reconnect in a time that is Alfvenic. 

Modeling the “Whole Sun Month” Corona with Improved Thermodynamics 

In the previous annual report we described our modeling efforts for the Whole Sun 
Month (WSM) time period (Aug. 10-Sep. 8, 1996). A paper on this work has been published 
(Linker et al. 1999), and is included in the Appendix. For our first 3D computation with our 
improved thermodynamic model, we recomputed the WSM case. The new calculation results in 
more realistic coronal temperatures, as shown in Figure 3. The improved thermodynamic 
description in our MHD model opens up the possibility of modeling disk emission, just as we 
have previously done for polarization brightness (pB). The more realistic temperature obtained 
from the solution can be used to predict the abundance of the coronal iron species and produce 
“simulated” EIT images. Figure 4a shows an FeXV 284A EIT image on August 27, 1996, 
while Figure 4b shows a simulated EIT image that was developed for this time period using the 
plasma parameters from the MHD computation. The simulated image reproduces the extension 
of the coronal hole past the solar equator, and also the small coronal hole in the south. The 
width of the dark emission region is much wider in the simulated image than in the EIT image; 
this may in part be due to the limited longitudinal resolution that was used for the computation. 

The active region also does not appear very bright in the simulated image, which may indicate 
that the coronal heating function used was too simple. Further comparisons of this type can 
help us to improve the MHD model as well as constrain models of coronal heating. 

Including the Transition Region in 2D MHD Models of the Solar Wind 

The simplified energy equation in polytropic MHD models fails to reproduce the 
temperature structure of the corona and the observed contrast in speed between the fast and 
slow solar wind. We have recently improved the energy equation in our MHD model by 
including thermal conduction parallel to the magnetic field, radiation, coronal heating, and 
Alfven wave pressure. We have also extended our model to include the transition region in our 
coronal calculations. Our lower boundary is placed at the top of the chromosphere (at a 
temperature of 20,000°K), so that the transition region is included in the domain of calculation. 

We specify a magnetic flux distribution on the solar surface and we integrate the time- 
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Figure 2. Reconnection of two current-carrying flux bundles. The time it takes 
to reconect the magnetic flux is on the order of the Alfven time. 
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3D MHD Calculation with an Improved Energy Equation for Whole Sun Month 

Temperature Field Lines Temperature 


Field Lines 


August 15, 1996 (CMD = 90 deg) 


August 22, 1996 (CMD = 0 deg) 


0.5 1.0 1 .5 2.0 3.0 MK 


Figure 3. The temperature in the corona for the model with an improved energy equation. Note that the 
closed-field regions (streamers) have higher temperatures, and the polar coronal holes have low temperatures. 

Simulated EIT Image of the Solar Corona 
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Figure 4. A comparison between (a) observed EUV emission at 284A from the SOHO EIT telescope and (b) a 
simulated emission image using plasma parameters from the MHD model. 

MHD Calculation of the Solar Wind with a Transition Region 


1 .5#4 K 2.0M K 


Figure 5. The plasma temperature (T) in the solar corona from an MHD simulation that includes the upper 
chromosphere and transition region. Also shown are magnetic field lines. Blue shows the lowest temperatures 
and red the highest. T varies from less than 20,000 K in the upper chromosphere to more than 2,000,000 K in 
the corona. 










SAIC-00/8005:APPAT-239 
March 22, 2000 

dependent MHD equations to steady state. The resulting solutions can be used to compare with 
the properties of the corona and of the solar wind. In this model, the solar wind is described 
self-consistently from its origins in the top of the chromosphere all the way out as it expands 
into interplanetary space. 

We have performed 2D (axisymmetric) MHD simulations of the solar corona, including 
the transition region, for a dipole magnetic flux distribution. This is a challenging computation 
because quantities vary by many orders of magnitude across short length scales in the transition 
region (down to 0.0005R 5 ). To model the transition region, we have incorporated several 
enhancements to our MHD code, including a more sophisticated treatment of the characteristic 
equations at the solar boundary, and an implicit advancement of the radiation loss term in the 
energy equation. 

Figure 5 shows the temperature at increasingly smaller scales with superimposed 
magnetic field lines. The temperature maximum of about 2 x 10 6 K is in the closed field 
regions under the streamer, where hot plasma is trapped by the magnetic field. Notice that the 
simulation is able to capture the sharp gradients in the transition region. A projection in 
Cartesian coordinates shows that the height at which the transition region forms varies with 
latitude. We presented a poster paper at the American Astronomical Society Meeting (Solar 
Physics Division), held in Chicago, Illinois, May 3 1-June 3, 1999, on “MHD Modeling of the 
Solar Wind Including the Transition Region,” by R. Lionello, J. A. Linker, and Z. Mikic. 

Modeling the Eruption of Arcades by Changes in Magnetic Flux 

Magnetic structures of various geometries, including loops and arcades, are present in the 
solar corona. Observations indicate that the magnetic field in some of these structures can be 
highly sheared, implying that a substantial amount of non-potential field energy is stored in the 
structure. If there is a physical mechanism that can induce a transition to a lower-energy state, 
the magnetic energy can be released into kinetic energy of plasma motions or thermal energy. 

We have studied the interactions between highly sheared structures (loops and arcades) and an 
emerging potential field structure by 3D numerical simulations. We found that the emerging 
field can destabilize an existing configuration, leading to the release of magnetic energy into 
plasma kinetic energy. A specific example is the eruption along the neutral line of a long, 
narrow, sheared arcade, which can be used to model a prominence eruption or a coronal mass 
ejection. 

We presented our progress at a contributed talk and a poster paper at the American 
Astronomical Society Meeting (Solar Physics Division), held in Chicago, Illinois, May 3 1-June 
3, 1999. The talk was on the “Initiation of Coronal Mass Ejections by Changes in 
Photospheric Flux,” by Z. Mikic and J. A. Linker. The poster paper was on “Eruption of 
Magnetic Structures in the Solar Corona,” by Y. Mok, Z. Mikic, and J. A. Linker. 
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Publication on Modeling of the Magnetic Structure of Prominences 

A paper that explores the three-dimensional magnetic field topology for prominence 
support has been published in the Astrophysical Journal (Amari et al. 1999), and is included in 
the Appendix. 

Presentation at the “Workshop on Physics of the Solar Corona and Transition 
Region,” Monterey, California 

We have investigated the effect of changes in photospheric magnetic fields on the stability 
of helmet streamers and active region arcades. Changes ih the magnetic flux in the vicinity of 
the neutral line can lead to disruption, with liberation of a significant fraction of the stored 
magnetic energy. When the amount of emerged flux is below a threshold, a stable equilibrium 
with a filament results. When the threshold is exceeded, the configuration erupts and leaves the 
Sun with a substantial amount of kinetic energy. This mechanism is a promising trigger tor 
launching CMEs. The results were presented in an invited talk at the Workshop on the Physics 
of the Solar Corona and Transition Region, held in Monterey, California, from 24—27 August, 
1999, in the talk “Modeling the Evolution of the Coronal Magnetic Field,'’ by Z. Mikic, J. A. 
Linker, and T. Amari. 

Reconstructing Force-Free Fields in the Corona 

During a visit by Dr. Tahar Amari (Ecole Poly technique, Paris) to SAIC in July-August 
1999, we collaborated on techniques to deduce coronal magnetic fields from vector 
magnetograph measurements of the magnetic field in the photosphere, using the force-free 
assumption. A paper on the mathematical aspects of a technique developed by Dr. Amari 
(Amari, Boulmezaoud, & Mikic 1999), with our collaboration, is included in the Appendix. 

Emergence and Interaction of Coronal Loops 

Recent observations from TRACE and SOHO indicate that magnetic loops are ubiquitous 
in the solar atmosphere. This discovery underscores the importance of understanding the 
dynamics and the potential consequences of these loops in the corona, as it is well known that 
dynamical events, such as solar flares, are often associated with these loops in active regions. 
Furthermore, how multiple loops form and coexist in a neighborhood is of fundamental 
importance because their interaction among each other could lead to the change of magnetic 
topology, resulting in the release of stored magnetic energy, which is believed to be the energy 
source of many dynamical phenomena. In our previous work, we have demonstrated a 
mechanism that leads to the dynamic formation of these loops, namely, by vortex plasma flows 
on the photosphere (Van Hoven et al. 1995). We have extended this loop formation study into 
the direct emergence of current-carrying loops from the photosphere. The physical mechanism 
and numerical method have been developed to show that a current-carrying loop with the 
observed properties can be dynamically formed in our proposed model (Mok et al. 1997). By 
modeling the time-dependent, normal components of the emerging loop’s magnetic field and 
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current density on the surface, we showed that the coronal field responds dynamically by 
forming a rising current carrying loop as if it emerges through the surface. 

A natural extension of this study is into an environment with multiple loops, such as those 
observed recently. In order to understand the dynamics of the multiple-loop configuration, we 
started by investigating a system of two loops, one loop emerging from the surface into a 
background with an existing overlying loop of comparable size and total flux in the 
neighborhood. This situation is not uncommon in an active region, where magnetic flux 
emerges frequently in the form of flux ropes. Using the same technique that we developed for 
the single-loop emergence, we simulated the formation of the second loop in a number ot 
configurations. Due to the fact that the parameter space is quite large even for a two-loop 
combination as shown below, we have studied only the cases that we believe will have a potential 
impact on the environment. In our study, each loop has a dominant toroidal magnetic-field 
component, and a poloidal component generated by the toroidal current. We found that the 
interaction of a two-loop system has four major critical parameters: 

(1) The angle between the two loop planes: we have studied three angles, i.e., 18 degrees, 

45 degrees and 80 degrees. 

(2) The relative direction of the toroidal magnetic fields: they can be “parallel” (Bj B 2 > 

0), or “anti-parallel” (B 1 B 2 < 0). 

(3) The relative direction of the toroidal current densities: they can be parallel (J1J2 > 0, 
or anti-parallel (J1J2 < 0). The combination of (2) and (3) determines the relative 
sign of magnetic helicity of the loops, i.e., whether they have the same sign or opposite 
sign. 

(4) The aspect ratios of the two loops and the relative size between their magnetic and 
current minor radii. 

The interaction in some of these combinations releases more magnetic energy than others 
on a shorter time scale, resembling a solar flare, although most of the combinations do not result 
in a violent interaction. The simulations start with a single magnetic loop, which has emerged 
and settled into an equilibrium state as discussed above. The results of one the most dynamical 
cases are described as follows. The second loop emerges underneath with a time scale that is 
long compared to the Alfven transition time along the loop in order to preserve the time scale 
ordering in the corona. One of the most dynamic cases is when the two loops are at 45 degrees 
with respect to each other. The toroidal magnetic field of the two are in opposite directions, i.e., 

Bi B2 < 0, while their toroidal current densities are in the same direction, i.e., J 1 ■ J2 > O' 1 ° other 
words, the loops have opposite sense of magnetic helicity. The (magnetic) aspect ratios of the 
loops are 3.3 and 6.6 respectively, while the radii of the current channel are half of the (minor) 
magnetic radii. The total magnetic fluxes of the two loops are the same. The time profile of the 
kinetic energy in the system is shown in Figure 6. The second loop begins to emerge at 
t = 300, when the first loop has reached equilibrium. After the decay of some initial 
disturbances, the kinetic energy rises on a time scale of 40 normalized Alfvdn times, compared 
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to the rise time of 200 of the emerging field. Talcing the nominal values of B = 200 Gauss, 
n = 5 x 10 8 /cm-\ major-radius = 10 9 cm, the Alfven transit time along the loop is ~ 10 seconds, 
and the time scale of the kinetic energy burst is ~ 21 seconds. Although the rise time of the 
emerging field in the simulations is somewhat faster than observations indicate, we believe that it 
is sufficiently slow to keep it separated from the kinetic energy rise time, so that the simulations 
do not consume excessive computer time. The total kinetic energy at its maximum is 2.4 x 10- 8 
erg, an appropriate amount for a small flare. The kinetic energy is contributed mostly by the 
horizontal velocity components. This is indicative of fie.Id line reconnection in the location 
where the rising apex of the emerging loop collides with the existing flux of the overlying loop. 

Since the toroidal components of B from the two loops have nearly opposite direction, the 
sharply bent reconnected field lines carry the plasma away from the reconnecting site in the 
classical FKR geometry. The interaction of the two structures takes place rather early, i.e., as 
soon as the apex of the emerging loop rises into the corona and is still at a distance from the 
major axis of the first loop. The evolution of the structure is shown in the sequence of 
Figure 7. The field lines gradually settle into a three-loop system. This interaction is a viable 
mechanism for the formation of multiple-loop configurations as recently observed in TRACE. 

The time profiles of several other cases are also shown in Figure 6. They typically evolve 
slowly into a new configuration with three to four loops connecting the four magnetic poles on 
the surface. 

Although a solar flare is a highly dynamical event and requires kinetic theory for a full 
treatment of the physical processes, such as particle acceleration, we have demonstrated that a 
substantial amount of stored magnetic energy is available to be released within the MHD 
framework. Exactly how the energy is converted into other forms, such as accelerated particles, 
X-rays, EUV, and microwaves, remains an open question and needs to be investigated. The 
results of this study are presently being prepared for publication (Mok, Mikic, & Linker 2000). 

Massively Parallel Version of the MAS Code 

During this year, we have spent a significant effort to port our spherical 3D MHD code 
(MAS) to massively parallel computers. This task was driven by the fact that high performance 
computing is being performed on massively parallel architectures. MAS is a FORTRAN 77 
code that was originally designed to am on vector supercomputers. 

The parallelization of a large scalar code can be accomplished in two fundamentally 
different ways: a) the code is entirely redesigned and rewritten from scratch, or b) the existing 
code is converted to the new architecture with the minimum number of modifications. After 
careful deliberation, we determined that the second approach would be most efficient. 

The Parallel Data Structure 

The parallel code was designed for a distributed-memory machine using the message- 
passing interface (MPI). The parallelization is obtained by considering a decomposition of the 
spatial domain and by assigning each sub-domain to a different processor. In other words, each 
processor stores in its memory and performs computations only on data that refer to a fraction 
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Figure 7a. Field lines at t = 325. The second loop (center) is emerging, and some of its field lines have 
reconnected to those of the outer loop. 



Figure 7b. Field lines at t - 350. More field lines from the lower loop have reconnected to those of the 
overlying loop, showing sharp turns in direction. 



Figure 7c. Field lines at / = 425. The field lines are settling into a three loop system. The original overlying 
loop is more or less gone. 
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of the overall domain considered in the simulation. The difficult part of the port was to “seam” 
properly the data on the different processors and to provide the required synchronization. 

As a general guideline, every variable in the old code becomes a “local” variable, in the 
sense that it refers only to the spatial sub-domain of each processor. The variables are still 
named in the same way but their value will be in general different on different processors. 
Naturally, a few exceptions are required and few new “global” variables, i.e., referring to the 
original whole domain, are introduced. These global variables are also defined on every 
processor, and their value is the same on every processor. • 

Implementing the Domain Decomposition 

The MAS code is essentially a time-dependent solver for a set of fields on a spatial 
domain, given specified boundary and initial conditions. This is a very general structure that 
applies to large variety of codes. There are two independent parts of the calculation: the time 
evolution and the field solution at each given time. The former is implemented with a time- 
stepping scheme, the latter by discretizing the field equations on a spatial grid. The 
parallelization concerns only the field solution, at each time step. Needless to say, nothing can 
be done to compute the time evolution in parallel; that is an intrinsically sequential calculation. 

The field equations in the MAS code are discretized with a finite-difference scheme on a 
structured mesh in spherical coordinates (/\0,0). A Fourier expansion is considered along the 0 
coordinate, leaving only a 2D mesh in the r-0 plane for the spatial domain decomposition. The 
number ot sub-domains along each of the two coordinates r and 6 is specified as an input; for 
example, by dividing the / -domain by 10 and the 0-domain by 6, the simulation will require 60 
processors (assuming that there is one sub-domain per processor). 

The domain decomposition must allow a certain overlapping between adjacent sub- 
domains to minimize the communication required between processors. The amount of 
overlapping is determined by the nature of the discretized differential operators that are used to 
represent the field equations. In the case of the MAS code, because a staggered mesh is used to 
discretize the equations, it turns out that a three-point overlapping along each linear dimension is 
required. For example, a domain decomposition into two equal sub-domains of a 40-point 
mesh along r will require the storage of mesh points 1-21 in the first sub-domain, and mesh 
points 19-40 in the second one. The three mesh points 19-21 are common to both sub- 
domains. 

Communication Among Processors 

At different times in the computation, data needs to be exchanged among processors. 

This communication has Been implemented by using MPI, a standard “message passing” 
library for parallel computing. All the global operations in the code require communication 
between processors. Here “global” refers to the whole simulation domain, as opposed to 
“local,” which refers to each sub-domain “owned” by each processor. Since communication 
time between processors is much slower than computation time on a processor, for optimum 
efficiency it is critical to reduce communication between processors to the minimum required. 
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The details of the MPI calls have been “hidden” in the code by writing the appropriate 
interface routines. Future improvements to the physics and/or to the numerical algorithms will 
be possible without requiring detailed knowledge of the inter-process communication 
procedures. 

Status of the Parallel MAS Code 

The parallelization of the MAS code has been fundamentally completed. The domain 
decomposition can be performed with an arbitrary number of processors along both the r and 9 
coordinates. Significant effort has been expended on the input/output interfaces to ensure that 
the input and output files are equivalent to the ones produced by the single-processor code. The 
parallel code has been validated by comparing with single-processor runs, and by varying the 
number of processors on multiple-processor runs. 

Presently, the efficiency of the scaling of the code on multiple processors is being 
investigated. Some accessory routines still need to be parallelized and/or revised to improve the 
parallel efficiency. We have found that the parallel scaling is strongly dependent on the number 
of mesh-points per processor: the greater this number, the smaller is the inter-process 
communication relative to the actual computation performed by each processor. More precisely, 
the “surface to volume” ratio of each sub-domain is in general the most significant figure to 
determine the impact of the data communication among processors. Optimization efforts, 
especially in regards to the communication routines, are presently in progress. We are 
presently still testing the code and will soon use it to do “production runs” with large meshes. 

Mapping in situ Solar Measurements Back to the Sun 

During this reporting period, our investigation has focused on mapping Ulysses and 
WIND in situ measurements back to the Sun using a combination of MHD models. We used 
a two-dimensional ( r,<p ), single-fluid, time-independent MHD model (Pizzo 1981) to map solar 
wind measurements back to the Sun as far as possible (typically 30-60 solar radii). This model 
makes the assumption that the measurements are in fact reversible, which may limit the 
applicability of the results. Test simulations, however, showed that the non-reversibility is 
localized in the vicinity of shocks and that overall, the solution may be considered “quasi- 
reversible.” The Pizzo model is applicable only outside of the outermost solar critical point. 
Furthermore, near the Sun there is considerable divergence of the flow in the meridional plane 
and thus modeling the flow in two dimensions becomes inappropriate. Therefore, to map the 
measurements from 30 solar radii to the solar surface we used our 3D solar coronal MHD 
model. Briefly, we traced field lines at the position of each of the mapped measurements (at 30 
solar radii) back to tbe solar’surface to deduce the source of the solar wind on the Sun. 

Figure 8 displays two views of these mapped measurements: (a) on the surface of a 
sphere; and (b) as a synoptic chart. In both cases, the poleward-most trace corresponds to 
Ulysses and the equatorward-most trace corresponds to WIND. The solar surface has been 
color-coded according to the topology of the magnetic field lines: black indicates open field 
lines (i.e., field lines connected to the Sun at one end only), and gray indicates closed field lines. 
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Figure 8. Ulysses and WIND trajectories are shown mapped back to the solar surface (for Carrington 
rotation 1913) and displayed (a) on the surface of a sphere; and (b) as a synoptic chart. The poleward-most 
trace corresponds to Ulysses and the equatorward-most trace corresponds to WIND. The solar surface has 
been colored according to the topology of the magnetic field lines; black indicates open field lines and gray 
indicates closed field lines. The trajectories are color-coded with their respective mapped speeds at 
approximately 30 solar radii; speeds 350 km/s and below are colored blue, speeds 700 km/s and above are 
colored red, and speeds in between are colored according to the color bar. 
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A striking feature ot this mapping is the significant variation in solar latitude, particularly 
at WIND. During this interval, Ulysses was located at a heliographic latitude of 27°N and 4.25 
AU, which maps back to a solar source latitude of approximately 70°N. In contrast, WIND was 
located in the ecliptic plane at l AU. These results, however, suggest that the source of the flow 
measured at WIND varied between 15°S and 70°N. WIND intercepted field lines from a small 
equatorial hole in the southern hemisphere and also intercepted flow from the equatorward 
extension of the northern polar coronal hole (the so-called '‘elephant’s trunk”). The mapped 
trajectories have also been color-coded with the inferred mapped speed at 30 solar radii. We 
remark that low speed (blue) maps back to the vicinity between the open (black) and closed 
(g ra y) field line regions. Moreover, for both spacecraft, speed tends to increase with distance 
away from this boundary. This is consistent with the view that fast solar wind originates in 
coronal holes and slow solar wind is associated with the boundary between open and closed 
field lines (Neugebauer et al. 1998; Linker et al. 1999). 

Modeling the Inner Heliosphere 

We have developed a three-dimensional MHD model to investigate the large-scale 
structure of the heliosphere (between 30 solar radii and 5 AU). Ultimately, we plan to drive this 
heliospheric model self-consistently using output from the coronal model. Currently, however, 
the simplified prescription of the polytropic energy equation employed in the algorithm does 
not yield sufficiently high plasma speeds necessary to drive a realistic heliospheric solution. 
Thus, as an interim solution, we utilize the magnetic field topology from the coronal solution to 
generate flow tields at the inner boundary of the heliospheric model. The coronal model itself 
is driven by the observed line-of-sight component of the photospheric magnetic field, and so the 
models can be, and are, run for specific time periods of interest. The results are compared with 
both remote solar observations as well as in situ observations and can be used as a basis for 
interpreting observations from a variety of disparate data sets. Figure 9 summarizes the 
heliospheric solution for the Whole Sun Month (a solar rotation involving parts of Carrington 
rotations 1912 and 1913). The heliospheric current sheet (inferred from the iso-surface B r = 0) 
is displayed out to 5 AU. The central sphere marks the inner boundary at 30 solar radii. A 
meridional slice of the radial velocity is shown at an arbitrary longitude. Blue corresponds to 
slowest speeds (350 km/s) and red corresponds to fastest speeds (750 km/s). Superimposed 
are a selection of interplanetary magnetic field lines, all emanating from the same heliographic 
longitude, but at different latitudes. Finally, the trajectories of the WIND and Ulysses 
spacecraft are displayed. A comparison of Ulysses and WIND measurements with the 
simulation indicates that the model has reproduced the essential large-scale features of the 
heliosphere during this time period. Figure 10 shows equatorial and meridional cuts of the 
solar wind velocity in the inner heliosphere during this time. 

Presentation at the Ulysses Science Working Team Meeting, San Diego 

In October 1999, Pete Riley presented a status report on our modeling efforts of the inner 
heliosphere, including the comparison with Ulysses in situ measurements. 
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The Heliosphere During Whole Sun Month 

August - September 1996 



Figure 9. Structure of the solar wind speed, magnetic field lines, and heliospheric current sheet in 
the inner heliosphere during the Whole Sun Month time period. 
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Figure 10. Structure of the solar wind speed in the inner heliosphere during the Whole Sun Month 
time period. Note the flow in the equatorial plane illustrates the characteristic "Parker spiral." 
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Presentation at the Western Pacific Geophysics Meeting, Taiwan 

In July 1998, Jon Linker presented an invited talk, “Magnetohydrodynamic Modeling of 
the Solar Corona For Space-Weather Applications,” at the Western Pacific Geophysics 
Meeting held in Taipei, Taiwan. 

Presentations at the Spring AGU Meeting and the SHINE 99 Workshop 

We presented the poster paper “Probing the Relationship Between the Solar Corona and 
Solar Wind Structure using MHD Models,” by P. Riley, Z. Mikic, J. A. Linker, J. T. Gosling, 
and V. J. Pizzo, at the Spring American Geophysical Union meeting, held in Boston, MA, May 
31-June 4, 1999, and at the SHINE 99 Summer Workshop, held in Boulder, Colorado, June 
14-17, 1999. We also presented the poster paper “Initiation of Coronal Mass Ejections by 
Changes in Photospheric Flux,” by Z Mikic, J. A. Linker, T. Amari, and J.-F. Luciani, and the 
poster paper “Relationship Between Electron Density and Temperature Within CMEs,” by P, 
Riley. J. A. Linker gave an invited review talk on “Existing Solar Wind Propagation Models” 
at the SHINE 99 meeting. 
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Magnetohydrodynamic modeling of the solar corona 
during Whole Sun Month 
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Abstract. The Whole Sun Month campaign (August 10 to Septembei 8, 1996) 
brought together a wide range of space-based and ground-based observations of the 
Sun and the interplanetary medium during solar minimum. The wealth of data 
collected provides a unique opportunity for testing coronal models. We develop a 
three-dimensional magnetohydrodynamic (MHD) model of the solar corona (from 
1 to 30 solar radii) applicable to the WSM time period, using measurements 
of the photospheric magnetic field as boundary conditions for the calculation. 

We compare results from the computation with daily and synoptic white-light and 
emission images obtained from ground-based observations and the SOHO spacecraft 
and with solar wind measurements from the Ulysses and WIND spacecraft. The 
results from the MHD computation show good overall agreement with coronal and 
interplanetary structures, including the position and shape of the streamer belt, 
coronal hole boundaries, and the heliospheric current sheet. From the model, we 
can infer the source locations of solar wind properties measured in interplanetary 
space. We find that the slow solar wind typically maps back to near the coronal 
hole boundary, while the fast solar wind maps to regions deeper within the coronal 
holes. Quantitative disagreements between the MHD model and observations for 
individual features observed during Whole Sun Month give insights into possible 
improvements to the model. 


1. Introduction 

The solar magnetic field plays a central role in coro- 
nal and interplanetary physics, defining the structure 
of the solar corona and inner heliosphere. Just as the 
number and complexity of sunspot groups vary with 
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the (approximate) 11-year solar cycle, the complexity 
of the Sun’s global magnetic field also varies, and this 
level of complexity is reflected in the structure of the 
solar corona. During solar minimum, solar activity is 
less frequent and the large-scale corona exhibits its sim- 
plest behavior. Long-lived helmet streamers and coro- 
nal holes may persist over several solar rotations. Solar 
minimum conditions provide us with an opportunity to 
separate the basic underlying structure of the corona 
from the solar active phenomena that disrupt it (such 
as coronal mass ejections), to understand how the so- 
lar magnetic field controls and influences the structures 
we see in the corona in various wavelengths, and to 
determine the solar source of phenomena measured in 
interplanetary space. 

With these goals in mind, a coordinated observing 
campaign of the solar minimum Sun was carried out for 
1 solar rotation. The campaign, known as Whole Sun 
Month (WSM), occurred from August 10 to September 
8, 1996, utilizing Solar and Heliospheric Observatory 
(SOHO), WIND, Ulysses, and ground-based data. The 
wide range of different solar observations collected al- 
lows us to obtain a more comprehensive understanding 
of the solar corona and its influence on the heliosphere 
during solar minimum. 

A fundamental difficulty in understanding coronal 
structure is that while the magnetic field is recognized 
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to play a crucial role, there are few measurements of 
the coronal magnetic field. However, the line-of-sight 
component of the magnetic field has been measured in 
the photosphere for many years (for example, at the 
Wilcox Solar Observatory and at the National Solar 
Observatory at Kitt Peak). A key to understanding the 
structure of the corona and solar wind is to accurately 
extrapolate the measured photospheric field outward. 
The simplest and most widely used techniques for ac- 
complishing this task are based on potential field mod- 
els (e.g., the potential field source-surface model [ Schat - 
ten et a/., 1969; Altschuler and Newkirk\ 1969] and the 
potential field current sheet model [ Schatten , 1971]). 
While these models do not directly describe the coronal 
plasma, they can address many aspects of coronal and 
interplanetary data [e.g., Hoeksema et a/., 1983; Wang 
and Sheeley , 1988, 1992, 1995], and the models are still 
being refined and improved [e.g., Zhao and Hoeksema , 
1995], Zhao et ai [this issue] discuss source-surface 
modeling for the WSM time period. 

Ideally, a model should not only extrapolate the mag- 
netic field, but also should self-consistently describe 
the plasma as well. Magnetostatic models [Bogdan 
and Low , 1986; Bagenal and Gibson , 1991; Gibson and 
Bagenal , 1995; Gibson et al., 1996] include the plasma 
force balance in the lower corona (neglecting flow terms) 
and have been used to quantitatively model specific 
time periods (with assumptions of azimuthal symme- 
try). Gibson et al. [1997] and Gibson et at. [The three- 
dimensional coronal magnetic field during Whole Sun 
Month, submitted to Astrophysical Journal , 1998] de- 
scribe magnetostatic modeling for WSM. 

Full magnetohydrodynamic (MHD) computations can, 
in principle, describe both the lower solar corona and 
solar wind. One advantage to this approach is that 
no assumptions are made about the nature of coro- 
nal currents or other properties of the solution; an ini- 
tial, value problem is specified and the solution is com- 
puted. While MHD computations have been performed 
since the 1970s [e.g., Endler, 1971; Pneuman and Kopp , 
1971], they have usually addressed idealized proper- 
ties of the corona and not specific observations. Re- 
cently, MHD computations that incorporate observed 
photospheric fields into the boundary conditions [Mikic 
and Linker , 1996; Usmanov, 1996] have been used to 
model specific time periods and have been compared 
with eclipse and interplanetary observations [Mikic and 
Linker, 1996; Linker et al. , 1996; Linker and Mikic , 
1997]. In principle, MHD modeling is a very powerful 
approach for studying the solar corona and solar wind 
when a coordinated set of observations is available, as a 
comprehensive model is then required to synthesize the 
different measurements into a coherent picture. How- 
ever, the use of this approach to address specific ob- 
servations is relatively new and is still being explored. 
The WSM campaign provides an extensive number of 
observations with which to compare, from which we can 
determine the strengths and weaknesses of the model. 


In this paper, we present MHD solutions of the solar 
corona out to 30 R s (solar radii) for the WSM time pe- 
riod and we compare the results with available coronal 
observations. Section 2 briefly describes the methodol- 
ogy. In section 3, we show some of the basic proper- 
ties of the MHD solution we obtained, and we compare 
the global structure predicted by the model with white- 
light and emission line images, especially the structure 
of the streamer belt and the location of coronal hole 
boundaries. We focus on a particular feature of in- 
terest, the extended coronal hole known as the “ele- 
phant’s trunk” (named because of its elongated shape) 
that was observed in late August 1996. In principle, 
MHD solutions can be extended beyond 30 R s far into 
interplanetary space (for example, MHD solutions of 
the solar wind from 30/?., out to 1 AU have been per- 
formed [ Pizzo , 1991]), but with the present model the 
extra computational expense is probably not warranted, 
given the simplicity of the energy equation (as discussed 
in sections 2 and 3, more sophisticated thermodynamic 
processes must be included to accurately describe fast 
and slow solar wind streams). Therefore, to compare 
the model with available solar wind data, we take the 
approach of mapping solar wind features back into the 
computational domain, as has been done previously 
with potential field models (see Neugebauer et al. [1998] 
for a comparison of several models). Using this tech- 
nique, we describe predictions by the model for the so- 
lar source of solar wind features observed at 1 AU by 
WIND and at 4 AU by Ulysses (1 astronomical unit 
(AU) = 1.49 x 10 s km = 214 solar radii). Section 4 dis- 
cusses our results in the context of present and future 
coronal modeling. 

2. Methodology 

To compute MHD solutions for the large-scale corona, 
we solve the following form of the equations in spherical 
coordinates: 


„ 47T 

YxB = — J 

(i) 

c 

1 OB 

( 2 ) 

= -VxE 

c at 

_ v x B 

( 3 ) 

E + = 7] J 
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dp 

t£ +v -('v)= 0 

( 4 ) 


fdv \ 1 

P ( + v • Yv J = -J x B - Yp+pg-f V(^pYv) (5) 

^ + V • (pv) = (7 — 1) (-pV v + S) (6) 

where B is the magnetic field intensity; J is the electric 
current density; E is the electric field; and v, p , and p 
are the plasma velocity, mass density, and pressure, re- 
spectively. The gravitational acceleration is g, 7 is the 
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Plate 1. (a) Composite synoptic radial magnetic field map developed from synoptic National Solar Observatory at Kitt Peak 
(NSO/Kitt Peak) maps for Carrington rotations 1912 and 1913 (see text). Blue (red) values indicate field directed into (out of) 
the Sun. (b) The radial magnetic field map used in the MHD calculation. This map was developed from the map in Plate la 
by interpolation and smoothing. The scaling in the two images is not identical but was chosen to emphasize the small-scale 
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Plate 2. Magnetic field lines from the MHD calculation superimposed on the radial velocity, corresponding to a central 
meridian longitude of 191° (September 3, 1996). The radial velocity has been color contoured in the plane perpendicular to 
the line of sight. Closed magnetic fields are shown in red, and open field lines are shown in black. The solar surface is shown 
in yellow. Red field lines that appear to be open actually close behind the plane depicting the velocity. 
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Plate 3. The radial velocity, number density, and temperature versus latitude and longitude from the MHD calculation for 
3 different heights in the corona, (a) r = 1.5/?,, (b) r = 2.5/?,, and (c) r = 10/?,. 
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ratio of specific heats, 77 is the resistivity, v is the viscos- 
ity, and 5 represents energy source terms. The method 
of solution of (l)-(6) has been described by Mikic and 
Linker [1994]; details of the algorithm are also described 
by Lionello et al. [1998] and Lionello et ai [Stability 
of algorithms for waves with large flows, submitted to 
Journal of Computational Physics , 1998]. In practice, 
we use the vector potential A, where B = Y x A, to 
implement the algorithm. Through the use of staggered 
meshes, the vector identity V - Y x A = 0 is preserved in 
discretized form, so that V*B = V- J = 0is satisfied 
to within round-off error throughout the computation. 

The calculation described here was performed on a 
85x81 x 64 nonuniform r, 0 , <f mesh, with A r & 0.013 R s 
near the base of the corona and A# « 1.6° near the 
equator. The <j> mesh was uniform, with A <j) = 5,6°. 
The simulation domain extended out to 30 R s . A uni- 
form resistivity 1 ] was used, corresponding to a resis- 
tive diffusion time tr = 4ttR 2 /t]C 2 = 400 hours. The 
Alfven travel time at the coronal base (ta = R 3 /Va) 
for |B| = 10 G (typical of the high-latitude fields) 
was 8.3 min (Alfven speed Va = 1400 km/s), so the 
Lundquist number tr/ta was ~ 2900. A uniform vis- 
cosity v was also used, corresponding to a viscous dif- 
fusion time (tv = R]/v) of 40 hours. A well-known 
problem with numerical MHD computations is that the 
achievable values of resistivity and viscosity are much 
larger than solar values [e.g., Mikic and Linker , 1994], 
Our two-dimensional calculations of helmet streamer 
equilibria [e.g., Linker and Mikic, 1995] indicate that 
solutions with larger diffusion typically have slightly 
larger closed field regions and smaller magnetic ener- 
gies than those with smaller diffusion, but the solutions 
are qualitatively similar. 

The term S in the energy equation (6) includes re- 
sistive and viscous dissipation (r/J 2 + vpVv : Yv). Be- 
cause we use enhanced values of the viscosity and re- 
sistivity relative to coronal values, the heating rate as- 
sociated with these terms may be unphysically large, 
and in the calculations presented here we set S = 0. 
Typically, we find that inclusion of resistive and vis- 
cous dissipation has little effect on the overall magnetic 
field solution [Mikic and Linker , 1994], With 5 = 0, 
(6) yields polytropic solutions; that is, d(pp~ 1 )/dt = 0 
(where d/dt is an advective derivative). (In this case, 
7 becomes the polytropic index and is not necessarily 
related to the ratio of specific heats.) These solutions 
have the advantage that relatively simple models can 
match many of the properties of the corona. However, 
values of 7 close to 1 (7 = 1.05 for the calculations 
reported here) are necessary to produce radial density 
and temperature profiles that are similar to coronal ob- 
servations; this reflects the fact that important thermo- 
dynamic processes have been omitted from the energy 
equation [Parker, 1963]. The inclusion of more sophis- 
ticated thermodynamics in (6), such as coronal heating, 
thermal conduction parallel to B, radiative losses, and 
Alfven wave dissipation, will be the subject of a future 


paper. These processes have been considered previously 
in one-dimensional (1-D) models [e.g. Withbroe, 1988; 
Habbal et ai, 1995], and some of these effects have been 
incorporated in two-dimensional (2-D) models [Suess et 
al, 1996; Mikic et ai , 1996a]. 

To perform a computation of the solar corona ap- 
plicable to the WSM time period, we obtained syn- 
optic maps (collected during a solar rotation by daily 
measurements of the line-of-sight magnetic field at cen- 
tral meridian) from the National Solar Observatory 
at Kitt Peak (NSO/Kitt Peak). Using the NSO/Kitt 
Peak maps for Carrington rotations 1912 and 1913, we 
developed a composite map using data from the pe- 
riod August 12 (CR1912, 122° Carrington longitude) 
to September 9, 1996 (CR1913, 122° longitude). This 
interval is centered on August 26 (CR1913, 303°), about 
the time that the “elephant’s trunk” coronal hole was 
at disc center. The composite NSO/Kitt Peak map was 
used to specify the radial magnetic field at the pho- 
tosphere B r 0 (in the manner described by Wang and 
Sheeley [1992]). Line-of-sight projection effects make 
accurate measurement of the polar fields difficult, so 
that the measured polar fields are less reliable. For the 
case shown here, we replaced the magnetic field within 
11° of each pole with values extrapolated from a smooth 
fit of the axisymmetric component of the measured field 
between 11° and 30° of the poles. A very small monopo- 
lar component present in the data was subtracted out. 
The NSO/Kitt Peak data are interpolated to the grid 
resolution used in the calculation, with care taken to 
preserve the magnetic flux in the interpolation. Plate 
1 shows the original NSO/Kitt Peak synoptic map, as 
well as the resulting map used for the MHD computa- 
tion. In the map used in the computation, fields near 
the poles had strengths of 10-12 G, while fields near the 
active region were ~ 40 G . 

In addition to B r0 , the MHD calculations require the 
specification of density p 0 and pressure p 0 at the coronal 
base [( Linker and Mikic [1997] describe the boundary 
conditions in more detail). We assumed AT, the plasma 
number density (that, is, N = n e = n;, where n e and m 
are the electron and ion number densities, respectively) 
at the coronal base to be 2 x 10 s cm’ 3 and the temper- 
ature To (defined by the ideal gas law p — 2 NKT) to be 
1.8 x 10 6 K. Rigid rotation corresponding to a sidereal 
period of 26 days was also imposed. 

For the initial condition for the fluid variables (p, 
p, and v), we choose the spherically symmetric tran- 
sonic Parker solar wind solution [Parker, 1963] consis- 
tent with the boundary values of p and p. In the ab- 
sence of a magnetic field and with no rotation, this so- 
lution yields a steady state. For the initial magnetic 
field, we compute a potential field by numerically solv- 
ing Laplace’s equation with boundary conditions on B r 
specified as B r 0 (shown in Plate 1) at r = R 3 and zero 
at the outer boundary. This yields an initial (current 
free) field of entirely closed field lines. Together, this 
specification of p , p , v, and B describes a nonequilib- 
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rium initial state for the time-dependent computation. 
We integrate the MHD equations (1-6) forward in time 
until a steady state is reached. The solution obtained 
provides a 3D description of the solar corona, under the 
assumption that the large-scale coronal magnetic field 
is not changing significantly for the time period of in- 
terest (an assumption that, while never strictly true, 
is often reasonably well satisfied during minimum solar 
activity). 

3. Results 

The typical solar minimum corona exhibits a well- 
defined streamer belt with coronal holes in the polar 
regions. The streamer belt is believed to be formed by 
closed magnetic field lines that have confined the solar 
plasma, while in coronal holes the solar wind streams 
outward freely along open magnetic field lines. In sec- 
tion 3.1, we show that the 3-D MHD solution we have 
obtained exhibits these basic properties. However, the 
primary purpose of our paper is to compare the solu- 
tions with available coronal and heliospheric data. This 
requires us to view the results in a manner that is similar 
to the actual observations. In sections 3.2 and 3.3, we 
describe comparisons of the MHD computations with 
white-light measurements and disc emission and absorp- 
tion images. The magnetic field structure predicted by 
the model can also be tested with interplanetary mea- 
surements. In section 3.4, we discuss comparisons of 
this type, and we use the MHD model to infer the source 
location of solar wind properties measured by WIND 
and Ulysses. 

3.1. MHD Solution 

We first show some of the genera] properties of the 3- 
D MHD solution. The properties of the plasma density 
and velocity from the solution are closely related to the 
associated magnetic field topology. This can be seen in 
Plate 2, which shows magnetic field lines from the MHD 
calculation superimposed on a color contour plot of the 
radial velocity v r in the plane perpendicular to the line 
of sight (v r is the dominant component of the velocity 
and closely resembles \v\). The plate shows that, as 
expected, the closed field regions are associated with 
nearly stagnant plasma flow. There is slow flow near 
the streamer belt in the equatorial regions, and faster 
flows are associated with the polar regions of the corona. 

The contrast of plasma properties between the open 
and closed field regions can also be seen in Plate 3, 
which shows u r , N (the plasma number density), and 
T (the plasma temperature) at r — 1.57?* (Plate 3a), 
r — 2.5 R s (Plate 3b), and r = 10 R s (Plate 3c). 
The streamer belt can be identified as the region of 
nearly zero flow velocity and increased plasma density 
at r — 1.5 R s and r = 2.5 R s (we demonstrate the closed 
magnetic field topology of this region in Figure 2). The 
strong “warp” in the structure of the streamer belt seen 
at 240° and 300° longitude is associated with northern 


and southern extensions of the coronal hole boundaries, 
which are discussed in subsequent sections. Higher N 
in the streamer belt at r = 1.5 R s is associated with the 
active region near 240° longitude. The plasma temper- 
ature also shows different properties between open and 
closed field regions; however, we caution the reader that 
because of the simplistic nature of the energy equation 
chosen, T is likely to be the least, realistic aspect of the 
calculation. Variations in T at each height are typically 
less than 10%, and isothermal ( T — constant) solutions 
would probably yield qualitatively similar properties to 
those shown here. 

At 10 R s , the magnetic field is essentially open every- 
where. Faster flow is still present at the poles relative 
to the equatorial regions, and the region near the helio- 
spheric current sheet is marked by higher A r than the 
surroundings. The low plasma pressure region (both 
low N and T near 240° longitude) is spatially coinci- 
dent with stronger magnetic fields arising in the ex- 
tension of the coronal hole here. At 10 R s and above, 
the flow is super-Alfvenic and supersonic nearly every- 
where. At this height and outward, the morphological 
features of the solution change very little, because in- 
formation propagates outward only. Quantitatively, N 
and T decrease with radius while v r increases, as is 
shown in Figure 1. 

The MHD solution shows the qualitative features we 
expect of the corona: stagnant flow and higher density 
in the closed field regions, faster flows near the poles, 
and slower flow near the streamer belt. To demonstrate 
that the solution can describe specific coronal observa- 
tions, we must cast the results in a form that is similar 
to the measurements. Two of the most common forms 
of coronal data are white-light and disc emission im- 
ages. Views of streamers on the limb and coronal holes 
on the solar disk together provide diagnostics of the 
three-dimensional structure of the corona that can be 
used to test coronal models. In the following two sec- 
tions we show how the MHD solution compares with 
these observations. 

Despite the qualitative similarities, the solution has 
important quantitative differences from the real corona 
and solar wind. In particular, the density contrast be- 
tween open and closed field regions is smaller than ob- 
served, and the fast solar wind speed in the simulation is 
too low compared to measurements. These inaccuracies 
in the solution are directly related to the simple energy 
equation assumed; improvements of these aspects of the 
solution are discussed further in section 4. 

3.2. Comparison With White-Light Images 

White light emitted from the solar corona is predom- 
inantly produced by electron scattering in the coronal 
plasma. The signal measured by a coronagraph is thus 
sensitive to the line-of-sight integrated electron density 
of coronal structures. When the limbs of the Sun are 
viewed in white light (with the photosphere occulted), 
the dense streamers show up as bright regions while 
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Figure 1. (a) Radial velocity, (b) number density, and 
(c) temperature versus radius at two different latitudes, 
for longitude 281° (corresponding to the west limb of 
Plate 2). The radial cut at 25° latitude passes through 
the closed-field region and then cuts across open field 
lines that bend toward the equator. The cut at 85° 
latitude is strictly in the open field region. 


the less dense coronal holes appear dark. Deducing the 
three-dimensional plasma density directly from corona- 
graph images is possible using tomographic techniques 
[Zidowitz et a/., 1996; Zidowitz , this issue]. On the other 
hand, it is relatively straightforward to use the plasma 
density from a model to develop a simulated coronal im- 


age. Frequently, the polarization brightness pB is the 
observed quantity. The distribution of pB in the plane 
of the sky is proportional to the line-of-sight integral 
of the product of the electron density and a scatter- 
ing function that varies along the line of sight [Billings, 
1966]. Plate 4 shows simulated pB images from the 
MHD model for 4 days during the WSM time period. 
Also shown are tracings of magnetic field lines for the 
same views. For comparison, pB observed with the 
Maun a Loa Solar Observatory’s Mark 3 K-coronameter 
(MLSO MKIII) are also shown. The streamer belt 
structure predicted by the model is qualitatively similar 
to the observations, although there are also significant 
differences. 

Another way of examining coronal structure during 
a solar rotation is to assemble the white-light measure- 
ments throughout the rotation at a particular radius 
into a chart of brightness versus latitude and longitude. 
These are often called synoptic images. An example for 
MLSO MKIII data taken during WSM at r = 1.15 R s 
and r — 1.75 R $ is shown in Plates 5a and 5b. Plates 
5c and 5d (same format as Plates 5a and 5b) show data 
at r = 2.35 R s and r = 5 R s obtained from the large- 
angle and spectrometric coronagraph (LASCO) aboard 
the SOHO spacecraft. Also shown in Plates 5a-5d are 
simulated synoptic images from the MHD model, de- 
veloped by processing simulated pB images in the same 
manner as the MLSO MKIII and LASCO data. 

The simulated and actual synoptic images show many 
similar overall features, indicating that the structure of 
the streamer belt during WSM has been captured by 
the MHD model. With increasing height, the streamer 
belt thins because of the expansion of the magnetic field 
from the coronal holes. At 5 R s , the streamer belt be- 
gins to have a sheet-like appearance similar to the he- 
liospheric current sheet. Both the observations and the 
simulation show a dark region breaking up the streamer 
belt near 270°; this phenomenon is associated with the 
elephant’s trunk coronal hole, discussed in more detail 
in the next section. 

3.3. Comparison With Disc Images 

Emission from the disc of the Sun provides another 
means of studying coronal structure. Coronal holes are 
typically identified from disc images of emission (at vari- 
ous wavelengths) by their low emission intensity relative 
to other regions of the solar atmosphere. Coronal holes 
are believed to be associated with open magnetic field 
regions. The expansion of the solar wind from coro- 
nal holes is thought to give rise to the lower plasma 
densities and temperatures responsible for the observed 
lower levels of emission. Absorption in the 10830 A (He 
I) line is enhanced (that is, brightness is decreased) by 
increased EUV and X ray emission, so the foot points 
of bright X ray loops appear dark and coronal holes 
appear bright in He 10830 A images [Harvey, 1996]. 
The National Solar Observatory at Kitt Peak publishes 
daily maps of coronal hole boundaries deduced from He 
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10830 A images and measurements of the photospheric 
magnetic field polarity. There is often a good (but not 
exact) correspondence between coronal holes identified 
in He 10830 A and soft X rays [ Harvey , 1996]. 

Perhaps the most conspicuous coronal feature ob- 
served during WSM was the elephant’s trunk coronal 
hole that extended from the northern polar regions to 
the equator. The elephant’s trunk was apparent in 
several different wavelengths, including SOHO extreme 
ultraviolet imaging telescope (EIT) images, NSO/Kitt 
Peak He 10830 maps, and Yokhoh soft X rays. This 
interesting feature was most visible around August 26- 
27 and provides another test of the model’s ability to 
predict coronal structure. 

To determine the boundary between open and closed 
field regions, we traced magnetic field lines from each 
pixel on the lower boundary. By delineating the regions 
where the field lines returned to the Sun (closed-field 
regions) from the regions where the field lines reached 
the outer boundary (open-field regions), we produced 
a map of the coronal hole boundaries predicted by the 
MHD model, under the assumption that open-field re- 
gions correspond to coronal holes observed in emission 
(or absorption) lines. Plate 6 show r s the results of this 
calculation. Plate 6a shows tracings of the magnetic 
field wfith B r mapped on the solar surface, as w r ould 
have been seen on August 27, 1996. Plate 6b show r s 
the same field tracings and the same view but with 
the open-field regions (black) and closed-field regions 
(gray) mapped on the surface of the Sun. For compar- 
ison, the NSO/Kitt Peak coronal hole boundaries and 
SOHO EIT images showung the elephant’s trunk coro- 
nal hole are shown in Plates 6c-6e. It is apparent that 
the mode] has captured the extended coronal hole fea- 
ture, although the observations show the coronal hole 
extending to even lower latitudes and hooking eastward 
into the active region. 

We note that the model also predicts an equatorial 
coronal hole in the south (east of the elephant’s trunk). 
A coronal hole in the south appears to be present in 
the EIT 195 A (Fe XII) and 284 A (Fe XV) images, 
though it is partially obscured by the bright emission 
from the nearby active region. However, it has not been 
identified as a coronal hole in the NSO/Kitt Peak map. 
This highlights one of the difficulties with trying to in- 
terpret magnetic field topology from spectroscopic mea- 
surements: a coronal hole can appear to be present in 
one wavelength but may not be visible in another. In 
the synoptic white-light image (Plate 5) a break appears 
in the streamer belt extending from both the north and 
the south. This feature is present in both the observa- 
tions and the MHD computation and, according to the 
model, is caused by the northern and southern coro- 
nal holes “pinching” on either side of the active region. 
This dark region in the white-light data suggests that 
the southern equatorial coronal hole is indeed present. 
As discussed in section 3.3, the detection of a polarity 
reversal in the interplanetary magnetic field that maps 


back to the southern coronal hole extension further sup- 
ports this interpretation. 

Plate 7a shows a “synoptic” EIT 195 A image, gen- 
erated by assembling measurements taken at central 
meridian, similar to the construction of synoptic charts 
of the magnetic field. The measurements have been 
binned by sine(latitude), which has the effect of com- 
pressing the polar regions into a smaller area. The im- 
age shows the polar coronal holes and elephant’s trunk 
for the WSM time period. Note that the coronal hole 
area in the northern polar regions is larger than in the 
south. This asymmetry in the coronal holes is enhanced 
because of the 7° tilt of the solar rotation axis present 
during this time period, but would still be visible even 
if there wore no tilt. The southern equatorial coronal 
hole (visible in the daily images, as discussed above) 
is more difficult to see in the synoptic image, perhaps 
because it is obscured by the neighboring active region. 

Plate 7b show r s the coronal hole boundaries predicted 
by the MHD calculation, in the same sine(latitude) for- 
mat as Plate 7a. Open-field regions are depicted as 
black, and closed-field regions are shown in gray. There 
are many qualitative similarities between the coronal 
holes identified in EIT and the open-field regions in the 
MHD model, including the presence of the elephant’s 
trunk coronal hole and the larger area of northern coro- 
nal holes relative to southern. However, the coronal 
hole area in both the north and south is larger in the 
MHD model than in the EIT image. The most notable 
disagreement between the model and observations are 
the equatorial coronal holes in the north near 40° and 
350° Carrington longitude. While a southward exten- 
sion of the coronal hole boundary seems to be present 
in the EIT data near 40° , it is not nearly as pronounced 
as in the MHD calculation. Neither the EIT daily im- 
ages nor the NSO/Kitt Peak daily coronal hole maps 
shows obvious equatorial coronal holes in these regions. 
These features are perhaps the most significant differ- 
ence between the MHD model and Whole Sun Month 
observations, and w^e discuss them further in section 4. 

We plot the open-/closed-field boundaries as a func- 
tion of radius in Figure 2. Note the close association 
between the field topology at 7’ — 1.5 R s and r = 2.5/tb 
and the plasma properties shown in Plate 3. At 5 R s 
the closed-field region resembles the heliospheric cur- 
rent sheet, and above this radius the field is essentially 
open everywhere. The magnetic field expands rapidly 
from the coronal holes; therefore, as showm in the next 
section, field lines that are located near the equator 
at large distance from the Sun frequently map back to 
much higher latitudes at the solar surface. 

3.4. Inferring the Source Location of In Situ 
Solar Wind Measurements 

One of the reasons for intense interest in coronal holes 
is because they are believed to be regions of open mag- 
netic field where the solar wind expands outward. The 
solar wind is far from uniform and is often described 
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3D MHD Model (B r ) 3D MHD Model (Coronal Holes) 



NSOKP Coronal Hole Map 



EIT Fe Xn Image EIT Fe XV Image 


Plate 6. Comparison of the MHD model with coronal holes deduced from measurements of emission and absorption lines, 
(a) Magnetic field lines from the MHD model, as would be viewed on August 27, 1996. B# is contoured on the solar surface 
in the same format as Plate 4c. (b) Same as Plate 6a, with open-field regions shown in black and closed-field regions shown 
as gray. Assuming that open-field regions in the model correspond to coronal holes, a feature similar to the elephant’s trunk 
coronal hole is produced by the model, (c) A coronal hole map for August 27, 1996, from NSO/Kitt Peak. The elephant’s 
trunk coronal hole is clearly visible. Note that it actually extends below the equator, (d) An extreme ultraviolet imaging 
telescope (EIT) 195 A (Fe XU) image for the same day, TTie elephant’s trunk coronal hole is visible and very similar to the 
Kit t Peak coronal hole map. Note that a possible coronal hole is visible in the south, east of the active region, (e) Same as 
Plate 6d, but for an EIT 284 A (Fe XV) image. 
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Plate 8. (a) Magnetic field lines in the MHD model, together with the heliospheric current sheet (HCS) and the Ulysses and 
WIND trajectories (mapped back to 15/Q. The field lines traced from Ulysses (green) and WIND (cyan) are shown for data 
corresponding to the Whole Sun Month period. The Ulysses and WIND trajectories are shown in the Sun’s rotating frame, (b) 
A close-up view of the field lines near the Sun. The model correctly predicts that Ulysses did not cross the HCS during Whole 
Sim Month but that WIND crossed briefly below the HCS and then returned to regions of positive polarity. 
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Carrington Longitude 


Figure 2. Open (black) and closed (gray) field line 
regions, as determined from the MHD computation, as 
a function of radius. 


as being either '"slow” (300-500 km/s) or “fast” (600- 
800 km/s). The solar origin of fast and slow solar wind 
streams is poorly understood. Measurements from the 
Ulysses spacecraft have demonstrated that uniformly 
fast solar wind is present at the poles of the Sun dur- 
ing solar minimum [Phillips et a/., 1995]. Spacecraft 
near the ecliptic plane, such as WIND, typically mea- 
sure predominantly slow solar wind interrupted by oc- 
casional high-speed streams. Equatorial coronal holes 
have long been associated with recurrent high-speed so- 
lar wind streams near Earth [Bohlin, 1977; Crooker and 
Oliver , 1994; Harvey , 1996]; therefore the identification 
of features like the elephant’s trunk coronal hole is re- 
garded as particularly important for deducing the solar 
source of solar wind measured near the ecliptic plane. 
During the WSM time period, WIND observed predom- 
inantly slow wind but intervals of fast wind were in- 
deed present. On the other hand, the Ulysses spacecraft 
(moving slowly southward after crossing the north pole 
of the heliosphere in April 1995) reached low enough 
latitudes (« 28°) that intervals of slow wind were mea- 
sured. In fact, the WSM period corresponded to the 
first encounter of slow solar wind measured by Ulysses 
since its foray into northern latitudes in early 1995 [Ri- 
ley et a/., this issue]. Thus both WIND and Ulysses saw 
transitions between fast and slow wind during WSM. 

To investigate the solar origin of the features observed 
by Ulysses and WIND, we located the source regions 
of the solar wind at the Sun using the following tech- 
nique. Since the MHD calculation w 7 as performed only 
to 30 R S) we first mapped the locations from the space- 
craft to 15 R s (to be well inside the domain of calcu- 
lation), using a ballistic approximation in which the 
change in longitude is computed from the time inter- 
val required for a plasma parcel to travel from 15K* to 
the spacecraft location, using the measured in situ so- 
lar wind velocity; we then used the MHD field lines to 
trace back from 15 R s to the solar surface. (Varying the 
mapping distance between 5 R s and 30 R s caused indi- 
vidual footpoints to change, but did not change overall 
trends.) Plate 8a shows tracings of field lines, includ- 
ing the Ulysses and WIND trajectories (mapped back 
to 15 R s ) and the heliospheric current sheet (HCS) (de- 
fined as the surface B r =0). Plate 8b shows a close-up 
of the field lines near the solar surface. Note that the 
WIND trajectory intersects the HCS during this time 
period, but the Ulysses trajectory does not. Therefore 
the MHD model predicts that two HCS crossings should 
have been observed by WIND during this time period, 
and no crossings should have been observed by Ulysses, 
which is consistent with the spacecraft magnetic field 
measurements, as described next. 

Plate 9a shows the magnetic field polarity (deter- 
mined from 1 hour spacecraft data, averaged over 6 
hours) mapped back to the Sun, as well as the coronal 
hole boundaries deduced from the MHD model. The 
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top trace shows Ulysses measurements; the lower trace 
shows WIND measurements. The polarity of the mea- 
sured magnetic field was determined by comparing the 
magnetic field direction angle with the Parker spi- 
ral angle 4> P (determined from the measured solar wind 
speed). Points for which was within 45° of 
were assigned positive polarity (shown in red); points 
for which was within 45° of <F P — 180° were assigned 
negative polarity (shown in blue); for all other points, 
the spiral direction did not obviously correspond to in- 
ward or outward fields (shown in white). Interplanetary 
magnetic field measurements usually have a high level of 
fluctuations because of wave activity in the solar wind 
and sometimes from the passage of transients. Individ- 
ual measurements may not lie along the spiral direction 
for either polarity, or may even have the opposite po- 
larity, without the presence of an HCS crossing. Except 
for discrepancies of this kind, with a perfect model of 
the heliospheric magnetic field, the red points should 
always map back to positive polarity regions of the Sun 
and blue points should always map back to negative 
polarity regions. 

For Ulysses observations (top trace in Plate 9a), we 
see that the predominant positive polarity for the mag- 
netic field is both predicted by the model and mea- 
sured at the spacecraft (red points mapping back to 
positive-polarity northern polar coronal holes). This 
is not surprising, as the spacecraft was still at high 
enough latitude to be well away from the HCS dur- 
ing this time period. For WIND, the model predicts 
that the spacecraft passed below the HCS (into nega- 
tive polarity fields) near 240° Carrington longitude (at 
the Sun) and back above the HCS (into positive po- 
larity fields) near 220° longitude. These are the points 
mapping back to the southern equatorial coronal hole. 
This first HCS crossing was predicted quite accurately 
by the model (blue points in the southern equatorial 
coronal hole). The passage of the spacecraft back to 
positive polarity appears to have occurred somewhat 
later than that predicted by the model, as evidenced 
by the clustering of blue points in the northern coronal 
holes between 170° and 220°. 

An analysis of the W IND magnetic field and electron 
pitch angle data indicates that it is difficult to accu- 
rately identify the HCS crossing back to positive polar- 
ity fields. Indeed, WIND measured extended intervals 
of mixed polarity, during this time, prior to the return 
to positive polarity. There were inconsistencies between 
the crossings identified from the magnetic field data and 
the electron pitch angle data. The sharpness of the 
crossing from positive to negative polarity (“inbound”) 
identified in the WIND data and the less clearly identi- 
fied return to positive polarity (“outbound”) is consis- 
tent with the shape of the HCS shown in Plate 8a: the 
HCS has a steeper gradient on the inbound crossing and 
a “glancing” encounter with the HCS on the outbound 
crossing. This may explain the apparent disagreement 


between the MHD model results and the magnetic field 
polarity measurements. 

The WIND observations of negative polarity fields 
suggest that the southern equatorial coronal hole iden- 
tified in section 3.2 is indeed present. This is clearly il- 
lustrated by the field lines in Plate 8, which shows that 
the magnetic field maps back to the southern equato- 
rial coronal hole (negative polarity) when WIND crosses 
below the HCS. 

The overall good agreement between the magnetic 
field polarities observed by Ulysses and WIND and 
those predicted by the MHD computations suggests 
that it is possible to use the model to deduce the source 
locations of features measured in interplanetary space. 
Plate 9b is in the same format as Plate 9a, except that 
measured solar wind speed for the WSM time period is 
now plotted, with red indicating the highest solar wind 
speeds and blue the lowest. The top (mostly red) trace 
is for Ulysses, which observed predominantly' fast wind 
but measured some slow wind that maps back to longi- 
tudes near 240° on the Sun. WIND observed primarily 
slow wind, but there were two intervals of fast wind; 
one interval maps back to the elephant’s trunk coronal 
hole and the other to near 120° longitude. 

Plate 10 presents two views of the same solar wind 
source location mapping as shown in Plate 9b. Plate 
10a shows the view from Earth on August 27; Plate 10b 
shows a polar view. A general trend one sees in these 
plates is that slow wind velocity usually maps back to 
regions close to the coronal hole boundaries; fast wind 
typically comes from deeper within coronal holes. Ad- 
ditional comparisons with solar wind speed determined 
from interplanet ary scintillation measurements are con- 
sistent with this interpretation [fLreen et a/., this issue]. 

4. Discussion 

The data comparisons described in the previous sec- 
tions show that the MHD model gives a good overall de- 
scription of the solar corona observed during Whole Sun 
Month. The approximate shape of the streamer belt 
and coronal hole boundaries^ including the elephant’s 
trunk coronal hole, are reproduced by the model. Two 
HCS crossings for the WIND spacecraft were predicted 
by r the model and, in fact, observed by the spacecraft. 

Of course, many details of the observations do not 
match the computational results. For example, the ele- 
phant’s trunk coronal hole extends to lower latitudes 
than predicted by the computation. The tilt and size 
of streamers observed on individual days sometimes dif- 
fer from those predicted, and there are many fine-scale 
features in the observations that are not present in the 
computations. Also, even during solar minimum, the 
Sun’s magnetic field is always evolving and a synop- 
tic magnetic field map can yield, at best, an approxi- 
mation of the state of the corona. For example, daily 
images show significant structural changes in the ele- 
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Figure 3. A comparison of coronal hole boundaries 
(open-field regions) (a) predicted by the MHD com- 
putation and source-surface models with source-surface 
heights of R ss = (b) 2.2R S) (c) 2.5 R s , and (d) 2.7 R s . 
Open-field regions are shown in black, and closed-field 
regions are shown in white. The MHD computation and 
the source-surface model with R ss = 2.2 R s yield similar 
results. The coronal hole boundaries recede as R ss is 
increased, but the elephant’s trunk coronal hole (near 
300° Carrington longitude) and the southern equato- 
rial coronal hole (near 240° Carrington longitude) are 
robust features. Similar variations in the MHD result 
are likely to occur if p 0 and/or T 0 are increased (see 
text). 


The need for an improved thermodynamic descrip- 
tion in an MHD model is particularly pertinent when 
investigating the nature of coronal holes as observed 
in emission lines. The comparisons we have shown for 
coronal hole boundaries and the open-field regions pre- 
dicted by the MHD model are, by definition, qualitative. 
The dark emission features identified as coronal holes 
clearly are strongly associated with magnetically open 
regions on the Sun, but the observations do not directly 
measure field topology. As we have seen for Whole Sun 
Month, images in different wavelengths can yield dif- 
ferent interpretations of the coronal hole boundary. To 
make a quantitative comparison between coronal mod- 
els and coronal hole observations requires a more so- 
phisticated model that can be used to reproduce the 
emission measurement itself. For example, the emission 
lines observed by EIT arise from the excitation of iron 
ions, a trace species in the coronal plasma. The iron 
population in the corona is especially sensitive to tem- 
perature. The polytropic energy equation used here is 
inadequate for describing the temperature of the corona 
to the desired accuracy. In the future, we hope that with 
a more sophisticated treatment of thermodynamics in 
the energy equation, we can model the emission line 
measurements and thus obtain a clearer understanding 
of the true correspondence between the dark emission 
features on the Sun that are identified as coronal holes 
and regions of open magnetic fields. 

Finally, we remark on inferences regarding solar wind 
sources. One of the fundamental questions in solar and 
heliospheric physics is the origin of the fast and slow 
solar wind. Wang and Sheeley [1990, 1994, 1995] have 
argued that the speed of the solar wind is primarily de- 
termined by the rate of magnetic flux-tube expansion, 
with slow wind occurring where the expansion factor is 
large. We have used the magnetic structure computed 
with the MHD model to infer the locations of solar wind 
properties measured by interplanetary spacecraft. We 
find that the slow wind typically maps back to near the 
boundaries of open-field regions, with fast wind coming 
from deeper within the open-field regions. Neugebauer 
et ai [1998] found a similar result for a comparison 
with Ulysses and WIND data during the Ulysses fast 
latitude scan (February-April, 1995). The flux-tube ex- 
pansion factor is usually larger near the boundaries of 
open-field regions, so our results are not inconsistent 
with the Wang and Sheeley picture. However, by con- 
struction, our model will map any solar wind feature 
back to somewhere in an open-field region. Another 
idea that has been proposed for explaining the fast and 
slow nature of the solar wind is that coronal holes are 
the source of the fast wind, and time-dependent pro- 
cesses temporarily open closed-field regions and supply 
the slow wind. If this idea is correct, then a steady-state 
description of the corona will not suffice for determin- 
ing solar wind source locations. Thus our results neither 
confirm nor rule out either of these interpretations. 
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Plate 9. Mapping of interplanetary measurements back to the Sun using the MHD model. Mapped points are plotted versus 
heliographic latitude and Carrington longitude at the Sun, with coronal holes (open-field regions) predicted by the MHD 
computation shown in black and closed-field regions shown in gray. The polarity of the photospheric magnetic field in the 
coronal holes is also shown, (a) The magnetic field polarity measured by Ulysses and WIND. Red (b\ue) points indicate where 
positive (negative) polarity fields were measured. Ulysses observed predominantly positive polarity fields and no HCS 
crossings, consistent with the predictions of the MHD model. WIND observed an HCS crossing to negative polarity the maps 
to near 240° Carrington longitude, consistent with the MHD model (blue points mapping to the southern equatorial coronal 
hole with negative polarity). The crossing back to positive polarity fields occurs later than predicted by the model, as 
discussed in the text, (b) The solar wind speed mapped back to the Sun. Red shows the fastest speeds, and blue shows the 
slowest. The interval of slow wind measured by Ulysses maps to near the coronal hole boundary (near 240° longitude). Two 
intervals of fast wind were observed at WIND; one maps to the elephant’s trunk coronal hole and the other to a southward 
extension of the northern coronal hole boundary near 120° longitude. 
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Earth View on August 27, 1996 
fa) 



Plate 10. Solar wind speeds measured at Ulysses and WIND and mapped back to the Sun using the MHD model, using the 
same data shown in Plate 9b. Red shows the fastest speeds, and blue shows the slowest, (a) View as seen from Earth on August 
27, 1996. The coronal holes (open-field regions) are shown in black. The highest-latitude trace shows Ulysses measurements, 
(b) View from above the solar north pole. Slow wind typically maps to near the boundaries of coronal holes, while fast wind 
typically maps to deeper within coronal holes. 
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A three-dimensional magnetohydrodynamic model ot the global solar corona is described. The 
model uses observed photospheric magnetic fields as a boundary condition. A version of the model 
with a poly tropic energy equation is used to interpret solar observations, including eclipse images 
of the corona, Ulysses spacecraft measurements of the interplanetary magnetic field, and coronal 
hole boundaries from Kitt Peak He 10830 A maps and extreme ultraviolet images from the Solar 
Heliospheric Observatory. Observed magnetic fields are used as a boundary condition to model the 
evolution of the solar corona during the period February 1997-March 1998. A model with an 
improved energy equation and Alfven waves that is better able to model the solar wind is also 
presented. © 1999 American Institute of Physics . [S1070-664X(99)94805-X] 


I. INTRODUCTION 

The sophistication of models of the solar corona has in- 
creased considerably since the idealized models of the 
1 980’s (e.g., see Mikic 1 and Low 2 for examples of early 
models). This has been brought about by a confluence of 
three key elements. First, the collection of high-resolution 
observations of the Sun, both in space and time, has grown 
tremendously. For example, consider the changes in our per- 
ception of the Sun brought about by Yohkoh 2 images of the 
x-ray Sun; Ulysses measurements of the polar solar wind; 4 
high-resolution white-light movies of solar granulation; 
high-resolution vector magnetographs of active regions; 5 So- 
lar Heliospheric Observatory (SOHO) high-resolution im- 
ages of the photospheric magnetic field, white-light corona- 
graph, and extreme ultraviolet (EUV) images of the corona, 
and energetic particle and solar wind composition measure- 
ments. The space-based observations are nicely comple- 
mented by an extensive archive of ground-based observa- 
tions [in particular, magnetic field measurements at Kitt Peak 
National Solar Observatory (NSOKP) and Wilcox Solar Ob- 
servatory (WSO); He 10 830 A observations of coronal hole 
boundaries at NSOKP; Mauna Loa coronagraph images; and 
interplanetary scintillation (IPS) measurements]. Second, the 
power and availability of supercomputers has made two- and 
three-dimensional modeling routine; the increase in comput- 
ing power that massive parallelism promises will extend the 
possibilities for realistic modeling even further. Third, the 
sophistication of the models themselves, both in their geo- 
metrical realism and the physics that has been included, has 
matured significantly. 

The application of magnetohydrodynamic (MHD) coro- 
nal models to these solar observations Tias begun to exploit 
this confluence of capabilities. It is now possible to make 
direct comparisons between observations and models of the 
solar corona, as illustrated below. The development of this 
modeling capability is especially timely, since the observa- 
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lions from present missions [including the recently launched 
transition-region and coronal explorer (TRACE) mission] 
and coming missions [Solar-B, set for launch in 2004, will 
provide high-resolution measurements of vector magnetic 
fields, x-ray, and EUV emission in active regions; the Solar- 
Terrestrial Relations Observatory (STEREO) will take coro- 
nal images from multiple viewpoints; and Solar Probe will 
explore the inner solar corona] have challenged our under- 
standing of the Sun. 

To fully exploit the available data it is necessary to apply 
sophisticated models that use observational data as inputs 
and that produce observable quantities as outputs. Through 
this interplay ot observations and theory we can improve our 
understanding of the Sun and heliosphere. In this paper we 
show examples of how large-scale models of the solar co- 
rona can be used to make detailed comparisons with obser- 
vations. 


II. POLYTROPIC MHD MODEL 

A self-consistent description of the large-scale solar co- 
rona requires the coupled interaction of magnetic, plasma, 
and solar gravity forces, including the effect of the solar 
wind. For simplicity we first describe a “polytropic model,” 
in which an adiabatic energy equation with a reduced poly- 
tropic index y (i.e., smaller than 5/3) is used. 6 This is a crude 
way ot modeling the complicated thermodynamics in the co- 
rona with a simple energy equation. The primary motivation 
lor using a reduced y is the fact that the temperature in the 
corona does not vary substantially (the limit y— * 1 corre- 
sponds to an isothermal plasma). A typical choice, used here, 
is 1.05. Detailed comparisons of our results with coronal 
observations indicate that while this model matches many 
features of the corona (as shown below), it is not accurate 
enough to quantitatively reproduce the properties of the co- 
rona and solar wind. In particular, this simple model fails to 
reproduce the fast (-800 km/s) and slow (-400 km/s) solar 
wind streams 4 that are measured at Earth, nor does it repro- 
duce the contrast in density and temperature that is observed 
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between streamers and coronal holes. An improved model 
that does not suffer from these limitations is presented in 
Sec. IV. 

Theoretical arguments indicate that magnetic reconnec- 
tion is crucial to describe the structure and dynamics of the 
solar corona/' 9 Yohkoh observations also present strong 
evidence for the importance of magnetic reconnection. 10 We 
therefore include the effect of plasma resistivity (with the 
important caveat that numerical models require the resistivity 
to be enhanced compared to coronal values 11 ). In the resis- 
tive MHD model, the coronal plasma is described by the 
following equations: 

47 T 

V x B= — J, (1) 

c 

1 r)B 

V x E= (2) 

c dt 

I 

Ei — vxB= 7/J, (?) 

c 

(ip 

~r~ + V -(pv) = 0. (4) 

dt 

( d\ \ 1 

p | — + v-Vvj = — JxB— V/>- V/? u . + pg* V*(ppVv), 

(5) 

■~ + V-(/jv) = (y— 1)(— /)Vv + 5), (6) 

fit 

where B is the magnetic field, J is the electric current den- 
sity, E is the electric field, p, v. p, and T are the plasma mass 
density, velocity, pressure, and temperature, and the wave 
pressure p xv represents the acceleration due to Alfven waves 
(see Sec. IV). The gravitational acceleration is g 
= -g 0 f/r 2 , rj is the resistivity, v is the kinematic viscosity, 
and 5 represents energy source terms. The plasma pressure is 
p = {n e + n p )kT, where n e and n p are the electron and proton 
densities; for a hydrogen plasma, n e = n p . The polytropic 
model is defined by setting 1.05 and 5 = 0 in Eq. (6), and 
zero Alfven wave pressure in Eq. (5), p vv . = 0. 

We have developed a three-dimensional code (3D) 1 1,1 2 
to solve the MHD equations (l)-(6) in spherical coordinates 
This code has been used extensively to model the 
2D and 3D corona, including the structure of helmet 
streamers, 13-15 coronal mass ejections, 1 1,56-19 and the long- 
term evolution of the solar corona and heliospheric current 
sheet (see Sec. III). Related methods have been developed by 
Usmanov 20,21 and Pisanko. 22 The following boundary condi- 
tions are used. 1417,19 The radial magnetic field B r is specified 
at the solar surface r = R s (e.g., from synoptic magnetic field 
observations at NSOKP and WSO, or from full-disk magne- 
tograms). This field may evolve in time (see Sec. III). The 
boundary conditions on the velocity are determined from the 
characteristic equations along B. The plasma pressure, as 
well as the plasma density in regions where the radial veloc- 
ity is positive, are specified at r = R s . (In the calculations 
presented here, the boundary values of p and p were chosen 
to be uniform.) Characteristic equations are also used at the 


upper radial boundary, which is placed beyond the sonic and 
Alfven points (typically at r = 30/? J1 although we have per- 
formed simulations that have included Earth’s orbit in the 
domain. 19 and beyond). 

Typical parameters for the quiet corona -3 yield a value 
for the Lundquist number of 5~ 10 13 . Since it is not possible 
to perform wel 1-resol ved numerical computations 11 at this 
large value of 5, we must content ourselves with calculations 
at lower, but still substantial, values of 5~ 10 3 - 10 4 . We also 
use finite viscosity v corresponding to a viscous diffusion 
time r,,~~ 10 2 -J0 3 r A , where r v -R.]fv , t a ~R s !v a is the 
Alfven transit time, and v A is the Alfven speed. At the lower 
values of 5 used in the simulations we expect that the static 
current sheets in the solution (e.g., the heliospheric current 
sheet) will form in approximately the correct location, but 
will be broader than those in the corona. 

Pneuman and Kopp 24 developed the first 2D model of 
helmet streamer equilibria by solving the steady-state MHD 
equations. Our approach, and that used in many other calcu- 
lations, is to integrate the time-dependent MHD equations to 
steady state. 25-29 The following initial condition is used. For 
the given B r distribution at r= R s , a potential magnetic field 
(VxB = 0) is calculated in the corona, and a transonic 
spherically symmetric wand solution 6 is used to specify p , p, 
and v. The MHD equations are integrated in time until the 
plasma and magnetic fields settle into equilibrium. The final 
state has closed magnetic field regions (helmet streamers), 
where the solar wind plasma is trapped, surrounded by open 
fields (coronal holes), where the solar wind flows freely 
along magnetic field lines, accelerating to supersonic speeds. 
This model has also been used to study dynamic events in 
the corona, including coronal mass ejections (CMEs); 1 1,16-59 
this aspect will not be addressed in the present paper. 

A. Comparison with eclipse images 

Total solar eclipses offer an excellent opportunity to ob- 
serve coronal streamers. The white-light polarized brightness 
of the corona that can be measured during an eclipse can be 
simulated from our MHD solution by integrating the electron 
density along the line of sight in the plane of the sky (con- 
volved with a scattering function 30 and filtered with a radi- 
ally graded filter to mimic the effect of instrument “vignett- 
ing”). Our first attempt to model the corona 13,14 was 
performed subsequent to the eclipse of 3 November 1994, 
observed in Chile. Since then, we have made three predic- 
tions , before the actual eclipse date, using magnetic field data 
from the previous solar rotation, for the eclipses of 25 Octo- 
ber 1995, 19 seen in Vietnam and India; 9 March 1997, seen 
in Siberia, China, and Mongolia; and 26 February 1998, seen 
in the Caribbean. These predictions were published on the 
World Wide Web prior to the eclipses. We used NSOKP and 
WSO synoptic magnetic field maps for the calculations. The 
comparison of simulated polarization brightness with eclipse 
images is shown in Fig. 1. Comparisons with Mauna Loa 
MK3 coronagraph observations on several days during the 
solar rotations surrounding the eclipses have confirmed that 
the basic large-scale three-dimensional structure of the 
streamer belt has been captured in the model. The agreement 
between the model and the eclipse images is quite good. 
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Eclipse Comparisons 

Field Lines Polarizaliou U rightness liclipse Image 



October 24, 1905 



February 26, 1998 


FIG. 1. Comparison of MHD computations of the solar corona with total solar eclipse observations. The 3 November 1994 eclipse was modeled subsequent 
to the eclipse; the other three calculations were predictions, performed in advance of the eclipse date, using magnetic field data from NSOKP and WSO. AH 
the images are shown with geocentric north up. except for the 1997 eclipse, which has solar north up. (The sources of the eclipse images are listed in the 
Acknowledgments section.) 


especially considering that three of the cases are predictions. 
Note, however, that these eclipses occurred close to solar 
minimum, when the large-scale structure of the Sun changes 
slowly between solar rotations. We are planning to predict 
the state of the corona during the forthcoming total solar 
eclipse in August 1999, which will be seen in central and 
eastern Europe, the Middle East, and western Asia. This pre- 
diction will be our most challenging yet, since this eclipse 
will occur close to solar maximum, when the structure of the 
corona is expected to be considerably more complicated, and 
changing more rapidly, than for the cases we have simulated 
previously. 

B. Comparison with Ulysses measurements 

The coronal magnetic field not only defines the structure 
of the solar corona, but the position of the heliospheric cur- 
rent sheet (HCS), and the regions of fast and slow solar wind 
as well. Understanding how the Sun influences the structure 
of the inner heliosphere requires an accurate mapping of the 


photospheric magnetic field into the corona and beyond. 
Source-surface models 31 ” 38 provide predictions of the struc- 
ture of the magnetic field in the corona and heliosphere. 
Source-surface models are relatively simple to apply, and 
have yielded important insights into the structure of the he- 
liosphere, but a number of aspects of the Ulysses data are not 
described well by these models. 39-41 In particular, the latitu- 
dinal profile of the radial magnetic field and the extent of the 
HCS predicted by source -surface models show significant 
discrepancies from Ulysses observations. During May-June 
1993 the Ulysses spacecraft, which was located at 30 °S lati- 
tude, ceased to observe sector boundary (i.e., HCS) 
crossings. 39 The “classic” Wilcox source-surface 
model 34,35 predicted that Ulysses would cross the helio- 
spheric current sheet, whereas the MHD simulation correctly 
predicted no crossing. 13,19 The radial magnetic field from the 
MHD computation shows little latitudinal variation, consis- 
tent with Ulysses observations, in contrast to the source- 
surface model result. This agreement with Ulysses data indi- 










Gray: B r > 0, Hlack: B r < 0 

FIG. 2. Heliospheric current sheet for Carrington rotation 1892 and the 
Ulysses trajectory {February- April. 1995). in the rotating frame of the Sun. 
The color of the trajectory indicates the polarity of the magnetic field, as 
measured by Ulysses. Two views of the HCS are shown, (a) from above the 
HCS, and (b) from below the HCS. The large-scale polarity of the magnetic 
field is consistent with that predicted by the MHD model (positive above the 
HCS, negative below); the short-scale differences may be due to Alfven 
waves and to structures that have not been included in the limited-resolution 
calculation. The MHD calculation was performed up to 400 solar radii. 
Ulysses was at a distance of — 1.4 astronomical units (AU) from the Sun at 
this time. 


cates that MHD computations may provide a better way of 
mapping phenomena in the solar wind back to their origin in 
the solar corona. 

The Ulysses fast latitude scan (February- April, 1995) 
was a time period during which the Ulysses spacecraft tra- 
versed rapidly from southern polar latitudes of the Sun to 
northern polar latitudes, and offers another opportunity to 
test the MHD model. 19 Figure 2 shows that the HCS struc- 
ture predicted by the MHD model generally matches Ulysses 
measurements of the magnetic field polarity. A more detailed 
analysis of the HCS crossings 19 shows consistency between 
the measurements and the model. We also used the MHD 
model (by tracing magnetic field lines, back to the Sun, in 
conjunction with ballistic mapping in the region outside the 
computation) to deduce the solar origin of plasma observed 
at Ulysses. The results suggest that the fast solar wind gen- 
erally comes from deeper within coronal holes than does the 
slow wind. 42 

C. Whole Sun Month comparison 


MHD modeling is particularly useful for studying the 
solar corona and solar wind when a coordinated set of obser- 


FIG. 3. Comparison of the MHD model with coronal holes seen in disk 
measurements on 27 August 1996. The “elephant’s trunk” coronal hole 
(extending from the north pole past the equator) can be seen in both the 
simulation and the data. The NSOKP coronal hole map is deduced from He 
10 830 A images. The SOHO/EIT images are in EUV wavelengths. 

vations is available, since it can help to synthesize different 
measurements into a coherent picture. The Whole Sun 
Month campaign (WSM; 10 August-8 September, 1996) 
brought together a wide range of space and ground-based 
observations during solar minimum. Our MHD model was 
used to interpret coronagraph and EUV images, 15 * 43 inter- 
planetary scintillation measurements of the solar wind 
speed, 44 and the structure of corotating interaction regions 
(CIRs) as deduced from energetic particle measurements. 45 

The “elephant’s trunk” coronal hole, an equatorial ex- 
tension of the northern polar coronal hole, named on account 
of its shape, was perhaps the most conspicuous coronal fea- 
ture observed during WSM, and was apparent in several dif- 
ferent wavelengths, including SOHO Extreme Ultraviolet 
Imaging Telescope (EIT) images, NSOKP He 10 830 A 
maps, and Yohkoh soft x rays. It was most visible around 
26-27 August 1996. Figure 3 shows tracings of the magnetic 
field from the MHD model as they would appear on 27 Au- 
gust 1996, with coronal holes (i.e., regions with open field 
lines, colored black) and closed- field regions (gray) mapped 
on the surface of the Sun. For comparison, NSOKP coronal 
hole boundaries and SOHO/EIT EUV images are also 
shown. It is apparent that the MHD model reproduces the 
elephant’s trunk coronal hole, 15 although the observations 
show that the coronal hole extends to lower latitudes than 
predicted by the model. 

We also investigated the solar origins of features ob- 
served by the Ulysses and Wind 46 spacecraft during WSM. 
The slow solar wind maps back close to coronal hole bound- 
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aries, while the fast wind typically comes from deeper within 
coronal holes, a pattern similar to that seen during the Ul- 
ysses fast latitude scan (Sec. II B). The model also predicts 
HCS crossings by Wind (but not by Ulysses) during the 
WSM time period. 15 Wind HCS crossings similar to those 
predicted were in fact observed. 

III. MODELING CORONAL EVOLUTION 

In the previous section we used our model to find steady- 
state coronal solutions for a given distribution of photo- 
spheric magnetic field. This approach is limited to the study 
of the long-time properties of the solar corona. In reality, 
even if we neglect large-scale eruptions like coronal mass 
ejections, the corona is changing continuously, even during 
times of solar minimum. This changing structure is driven by 
changes in the photospheric magnetic field: active regions 
emerge and disperse continuously during the solar cycle. We 
have extended our model to incorporate the evolution of the 
photospheric magnetic field, so that we can now follow the 
evolution of the corona. This gives us the capability to study 
the long-term evolution of the corona (as shown below), the 
detailed evolution during a time period of interest (e.g., dur- 
ing Whole Sun Month), as well as the ability to study theo- 
retically the coronal consequences of changes in photo- 
spheric magnetic flux. 

When we seek steady-state solutions of Eqs. (l)-(6), we 
set the tangential component of the electric field at the 
boundary, E, o , to zero. This keeps ( B r at r = R s ) fixed in 

time. To make the flux evolve to match observed changes, it 
is necessary to specify a nonzero E, q . In general, E, q can be 
expressed as V,x^r + V,<f\ where ¥ and <I> are arbitrary 
functions (of 6 and <f>) and V, indicates tangential derivatives 
(in the 6-cf) plane at r~R s ). The potential changes E, o 
without changing the flux B r0 , and can be used to control the 
transverse magnetic field (i.e., the shear and the normal elec- 
tric current), whereas the potential 'T changes the flux. Since 
line-of-sight magnetograms do not provide information 
about the transverse component of the magnetic field, we 
only consider the case d> = 0 here. Note that a nonzero can 
be used to introduce shear into the field 18 and to match vec- 
tor magnetic field measurements. 47 The potential T' is ob- 
tained by solving the equation cV;^ = dB r !dt. Therefore, 
V is evaluated as new solar magnetic field measurements 
become available, specifying the evolution of E f(}t which is 
used as a boundary condition for the MHD equations. 

Thus, rather than computing a sequence of steady-state 
solutions for each set of magnetic field boundary values, our 
time-dependent MHD model now' represents the actual state 
of the corona corresponding to the evolving magnetic field 
measured on the surface of the Sun. We have used a se- 
quence of synoptic Kitt Peak magnetic field observations to 
study the evolution of the corona during the period 1 Febru- 
ary 1997-18 March 1998 (14 Carrington rotations). This 
time interval covers the beginning of the new' solar cycle, as 
the Sun emerges from solar minimum, and includes the 
emergence of high-latitude active regions. To model the evo- 
lution over a time interval of over a year is computationally 


prohibitive at present. In order to study the quasistatic evo- 
lution of the corona, we changed the photospheric magnetic 
field at a rate that w'as enhanced by approximately ten times 
compared to real time. 48 This approximation makes it impos- 
sible to study the detailed evolution of individual events, 
though it is still meaningful to study the quasistatic evolution 
of the large-scale structure of the corona. Figure 4 shows the 
evolution of the streamer structure, the coronal hole bound- 
aries, and the heliospheric current sheet during this time pe- 
riod. Note thejncrease in complexity of the coronal magnetic 
field as the 3un emerges from solar minimum. The output of 
this model can easily be compared with coronal observa- 
tions, as was demonstrated in Sec. II. 


IV. IMPROVED MHD MODEL 

Detailed comparisons of our model with observations, 
such as those presented in Sec. II, have forced us to confront 
the limitations of the polytropic MHD model. While the 
polytropic model matches many features of the corona, it 
does not reproduce the properties of the fast and slow solar 
wind or the large contrast in density and temperature be- 
tween streamers and coronal holes. We have improved the 
model by including the physical mechanisms that describe 
the transport of energy in the corona and solar wind. One- 
dimensional models have been quite successful, despite their 
obvious geometrical limitations, in describing this 
interaction. 49,50 We have improved the energy equation in 
our model to include the effects of parallel thermal conduc- 
tion, radiation loss, and parameterized coronal heating, and 
we have included a self-consistent model for Alfven wave 
acceleration. The source term in the energy equation (6) is 
given by 

S=-V-q-n v n p Q(T) + H* i + H d +D, (7) 

where // ch is the coronal heating source, D is the Alfven 
wave dissipation term, 77J 2 + rVvTv represents heat- 
ing due to viscous and resistive dissipation, and Q(T) is the 
radiation loss function. 51 In the collisional regime (below 
~ 10 R x ), the heat flux is q= - where 6 is the unit 

vector along B, and « , | j = 9 X 10 7 7 5/2 is the Spitzer value of 
the parallel thermal conductivity. The polytropic index y 
now becomes the ratio of specific heats, y= 5/3. In the col- 
lisionless regime (beyond ~~ 10/?,), the heat flux is modeled 
by q-an^Ty, where a is a parameter. 52 Since it is pres- 
ently not known in detail what heats the solar corona, 9 the 
coronal heating source /7 ch is a parameterized function. A 
typical form is 

H ch =H {) (6)cxp[-^-R s )/\(d)l (8) 

where H^(6) expresses the latitudinal variation of the volu- 
metric heating, and \{6) expresses the latitudinal variation of 
the scale length. [In practice, the variation can be expressed 
in terms of the magnetic topology (i.e., a proxy for the open 
and closed field regions) rather than the latitude 0.] 

Since the acceleration of the solar wind by Alfven waves 
occurs on spatial and time scales that are below the resolu- 
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Evolution of the Solar Corona During Feb. 1997 - Mar. 1998 


Coronal Holes 



CR 1919 



CR 1922 



CR 1925 



CR 1928 



CR 1934 


Photospheric B r Polarization Heliospheric 

and Field Lines Brightness Current Sheet 



RG. 4. The changing structure of the solar corona during the period February 1997-March 1998. Carrington rotations (CR) 1919-1934, as illustrated by 
coronal hole maps (longitude vs latitude, with gray/black indicating closed/open field regions), field line traces with the radial magnetic field shown on the 
surface of the Sun. polarization brightness, and the shape of the heliospheric current sheet. The HCS is shown up to 30 solar radii. The photospheric magnetic 
field was set as a time-dependent boundary condition on the 3D MHD simulation using NSOKP synoptic maps. 


tion of our global numerical model, the wave pressure p w is 
evolved using a WKB approximation 53 for the time-space 
averaged Alfven wave energy density e: 

cl€ 

— + V -F= D, (9) 

where F - ( 3/2v + ) e is the Alfven wave energy (lux, v A 
= B / \'4 7 rp is the Alfven speed, and p w =ef2. The Alfven 
wave velocity is \ A =±bv A ; in a multidimensional imple- 
mentation, it is necessary to transport two Alfven wave fields 
(waves parallel and antiparallel to B), which are combined to 
give €. The Alfven wave energy density € is related to the 
space-time average of the fluctuating component of the mag- 
netic field SB by €=(SB 2 )/4it. The dissipation term D ex- 
presses the nonlinear dissipation of Alfven waves in inter- 
planetary space and is modeled phenomenologically. 52 


The boundary r~R s is now chosen to be at the top of the 
transition region, at a given temperature (say T 0 
= 500 000 K). The density at r-R x is determined by balanc- 
ing radiation loss, thermal conduction, and heating within the 
chromosphere and transition region. 49 In this formulation, 
the only boundary conditions required from observations are 
on the radial magnetic field. Instead of specifying a nonuni- 
form temperature at the coronal base to express the observed 
variation of temperature between streamers and coronal 
holes, as would be required in the polytropic model, we now 
specify the distribution of coronal heating. By investigating 
how MHD solutions compare with observations it will be 
possible to test different coronal heating models, and, even- 
tually, when the coronal heating process is better understood, 
to relate the heating source to physical quantities. 

This formulation has been applied to a 2D (axisymmet- 
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ric) model of the corona. Extensive tests show that a nonuni- 
form heating profile, together with Alfven waves, can repro- 
duce the speed, mass flux, density, and temperature of the 
fast and slow solar wind at Earth. The coronal density con- 
trast is much improved compared to the polytropic model: 
coronal streamers are 5-10 times denser than coronal holes. 
The model is presently being extended to 3D. 

V. CONCLUSION 

The last decade has seen a marked increase in the so- 
phistication of models of the solar corona. Present models 
have improved geometrical realism, improved physics, are 
able to use solar observations directly as boundary condi- 
tions on the calculations (e.g., measured photospheric mag- 
netic field), and can model the evolving solar corona. We 
have presented comparisons between a 3D MHD model and 
observations of coronal and heliospheric structure. The 
agreement between simulated coronal structure and observa- 
tions is encouraging, implying that the models are mature 
enough for detailed analysis. Such 3D models will undoubt- 
edly find increased use in interpreting solar observations and 
developing up-to-date models of the solar corona. Yet it 
must be recognized that we have only begun to “scratch the 
surface” of what is possible. Despite these advances, some 
of the fundamental theoretical questions in solar physics re- 
main unsolved [e.g., what initiates CMEs and flares; what is 
the relationship between CMEs and flares; how do coronal 
magnetic structures emerge, evolve, and erupt; what heats 
the solar corona?]. Continued comparisons of model predic- 
tions and observations will help to answer these key ques- 
tions and provide new insights into the physics of the corona. 
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ABSTRACT 

The search for a background magnetic configuration favorable for prominence support has been given a great 
deal of attention for several decades. The most recent theoretical studies seem to agree that a promising candidate 
for the support of the dense and cooler prominence material, which fulfills several of the theoretical and obser- 
vational requirements such as twist, shear along the neutral line, and dips, is a magnetic flux rope. The most 
convincing models take an infinitely long periodic configuration that consists of a linear constant-a force-free 
magnetic field. These models, however, assume values of a that are close to its maximum possible value. In this 
Letter, we report our recent results, which show that it is indeed possible to produce a configuration that consists 
of a twisted magnetic flux tube embedded in an overlaying, almost potential, arcade such that high electric currents 
(and therefore values of a) are confined to the inner twisted magnetic flux rope. We present two MHD 
processes — corresponding to two different types of boundary conditions — that produce such a configuration. Our 
results show that the process associated variations of B z at the photospheric level by applying an electric field 
involving diffusion is much more efficient for creating a structure with more twist and dips. 

Subject headings: MHD — stars: coronae — stars: flare — stars: magnetic fields 


1. INTRODUCTION 

Prominences are sheets of cold material (about hundreds of 
times cooler and denser than the surrounding corona) supported 
by magnetic forces. They may either be quiescent structures 
outside active regions or participants in solar eruptive phenom- 
ena, such as coronal mass ejection or flares. Prominences are 
seen to lie above a neutral line separating regions of opposite 
polarity (Priest 1982). A global self-consistent model represents 
a complicated task because of the strong coupling between 
mechanical equilibrium and energy transport. Therefore, only 
two uncoupled problems have generally been considered so far, 
In the first class of problems the magnetic field is fixed while 
the energy equation is solved, while in the second one the global 
magnetic background configuration is sought for a given sim- 
plified energy equation or even without any energy equation 
(looking for a configuration favorable for prominence support, 
assuming that the prominence material just represents a small 
perturbation to the configuration). 

Magnetic flux ropes represent good candidates for the second 
class of problems and have therefore been sought. They may 
contain twist, shear along the neutral line*(Heyvaerts & Hag- 
yard 1991), and dips (locations of the magnetic configuration 
at which the concavity is favorable to prominence material 
accumulation). They also represent a natural way to account 
for the total solar (internal plus external) magnetic helicity 
content, which must be globally conserved and is found in the 
corona if helicity is injected from the interior of the Sun (con- 

' Also at Centre National de la Recherche Scientifique, Observatoire de 
Paris, Laboratoire de Physique Solaire et de V Heliosphere, F-92195 Meudon 
Principal Cedex, France. 


vection zone) into the corona (Low 1994). They must also play 
an important role by participating in helicity ejection (or re- 
distribution at “infinity” — in the solar wind) during eruptive 
events, via instability processes in which the flux rope is either 
passive (if the overlaying structure is disrupted) or active (if 
the flux rope is at the origin of a kink-type instability related 
to the amount of twist or electric current profile) (Amari, 
Luciani, & Mikic 1999). Indeed, prominences have been ob- 
served in several eruptive active regions (Hundhausen 1988, 
1994; Lekaet ai. 1993). 

For configurations invariant by translation along the x-axis, 
if one assumes the background magnetic field to be potential 
and the prominence to be modeled by a thin cold massive 
current sheet, then for various possible boundary value prob- 
lems (BVPs), no magnetic island can exist in the configuration 
(Aly, Amari, & Colombi 1990; Ridgway, Priest, & Amari 
1991a, 1991b), which excludes any flux tube-like topology. It 
was even possible to prove rigorously that with the same 
assumption of axisymmetry and symmetry against the vertical 
z-axis, no force-free configuration (linear or even nonlinear) 
that corresponds to a strictly bipolar photospheric boundary 
condition having a dip in the configuration can be accessible 
by slow quasi-static shearing motions (Amari et al. 1991), a 
result that actually excludes not only flux ropes or configura- 
tions having magnetic islands, but extends to a more general 
class. Relaxing the assumption of strictly bipolar boundary 
conditions in the case of a constant-^ linear force- free magnetic 
field defined in a half-strip, it was eventually possible (Amari 
& Aly 1992) to build a whole class of magnetic configurations 
having a flux rope topology for either bipolar (but allowing 
parasite boundary polarities), quadrupolar, or hexapolar con- 
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figurations. However, one major drawback of these models is 
that all these configurations are obtained for large values of a, 
close to the maximum permitted value (a„ ~ ttIL if L is the 
width of the strip). Aulanier & Demoulin (1998) substitute 
periodicity for axisymmetry (jr-invariant field) in our previous 
models to get configurations matching observed structures, but 
keeping the same strong limitation on the large values of a, 
which is unrealistically large even in the far field region. 

In this Letter, we report some new results that we have 
obtained in revisiting this problem of building a twisted flux 
tube in equilibrium in an overlaying arcade as the global back- 
ground magnetic configuration (above a bipolar photospheric 
region) for prominence support, using a three-dimensional 
MHD approach corresponding to tw o different types of bound- 
ary conditions. In the first one, purely ideal photospheric mo- 
tions are applied. In the second one, we apply an electric field 
that corresponds to diffusion of B : . This has been used to model 
the dispersion of active region flux into the network (Wang, 
Sheeley, & Nash 1991), and it actually appears quite reasonable 
to consider that in a long-term evolution of an active region, 
photospheric diffusion does occur (van Ballegooijen 1999). It 
is, however, worth noticing that none of these assumptions is 
indeed necessary for the rest of this Letter; these processes may 
be simply regarded as tools for building equilibrium configu- 
rations as the solution of MHD equations, in the same way 
that one could derive the equilibrium solutions in the previous 
models on pure analytical grounds. 

The Letter is organized as follows. In § 2, we present the 
three-dimensional MHD evolution (buildup of shear along the 
neutral line) of an initially arcade-like potential magnetic con- 
figuration. In § 3, we consider an ideal MHD process that adds 
up twist in a second phase to build a first type of twisted flux 
tube configuration, while in § 4 we present a diffusive MHD 
process that allows us to produce flux ropes with higher twist. 


points on the bottom of the box such that the normal component 
S.(jc, >\ 0) of the field on the boundary is preserved, which, as 
in Amari et al. (1996a, hereafter ALAT), corresponds to two 
large-scale parallel vortices rotating in the same direction with 
a maximum velocity not greater than r n = 10~ 2 (then v 0 is small 
compared to the Alfv£n speed i\ = I ). On the other faces of 
the box, we take the homogeneous Dirichlet condition v = 0, 
which is natural for a viscous plasma in contact with a non- 
moving wall. This type of velocity field corresponds to high 
shearing motions near the neutral line (very similar to those 
considered in Amari et al. 1996b in the axisymmetric case), 
while twisting motions are introduced as one move away from 
the neutral line on [z = 0}. 

The system is then evolved by solving the MHD equations 
that write in a simplified nondimensionalized form (see Amari, 
Luciani, & Joly 1999; ALAT): 


p — — —p(v 9 Vr) 4- (V x B) x B 
dt 

— Vp + V • (vpVv) + pg, 

~ = V x(i?xfl) — Y x (777), 

at 

dp _ 

— = -V(pt>), 

at 

— = — (r • V)p — I>(V * r) + H , 

r)/ 

j = V x B, 

V • B = 0, 


( 1 ) 

( 2 ) 

(3) 

(4) 

(5) 

(6) 


2. THE BUILDUP OF SHEAR 

In what follows, the corona is represented by the half-space 
D = U>0}, which is filled up with a slightly resistive and 
viscous plasma. For the numerical computations, D is approx- 
imated by the finite box {0 <x<L x , 0 <>’ < L v , 0 < z < Z..}, 
whose size is large compared to the characteristic spatial scale 
of the system (which is our reference length). Typically we 
choose L x = L y = L z - 40-60 and provide the domain with a 
nonuniform mesh (111 x 101 x 70 nodes). We first start at 
t — 0 with an initial potential (i.e., current-free) arcade-like 
configuration such that B[(x y \\ c); 0] = V^(j:,y, z), where 'k 
is the unique solution of a Dirichlet-Neuman Laplace BVP. Its 
normal derivative on {z = 0} corresponds to the distribution 
of B z . It is taken as two elliptic Gaussian distributions sym- 
metrically placed across the jr-axis; their centers are (0, ±0.8), 
their width is (&c = I, 6y = 2), and the normalization is such 
that the maximum value taken for B z is 1 . Oh the other (lateral) 
boundaries, we actually impose for ^ the value taken by the 
scalar potential corresponding to the unique current-free mag- 
netic field in the whole upper half-space {z > 0) whose distri- 
bution of B : on [z = 0} coincides with the normal derivative 
of ^ on the finite computational box and tends to zero at 
infinity. With this large box and this type of boundary con- 
ditions, the initial condition in the finite computational domain 
D mimics an arcade-like solution in the upper half-space 
fe>0). 

For t > 0, we impose a twisting velocity field to the foot- 


where B is the magnetic field and t\ p, and v are, respectively, 
the velocity, mass density, and kinematic viscosity of the 
plasma, while rj is the plasma resistivity. Since in coronal sit- 
uations 0 is small (of the order of 10“ 3 or even smaller), typical 
simulations were done with these values or with 0 = 0. As in 
.ALAT, in the latter case we are thus constrained to fix arbitrarily 
a mass density profile and of course to neglect the gravity term 
in equation (1). Here, we choose p = B 2 , which ensures a 
constant Alfvdn velocity. We have also tried other different 
density profiles (exhibiting a slower decrease with distance), 
but this does not lead to any difference in the results. Small 
values are used for the dissipation coefficients: v = 10~ 2 to 
10~ 3 for the kinematic viscosity and 17 = 0 for the first sequence 
of runs. The calculations of the evolution of the field are per- 
formed with the code METEOSOL, which uses a implicit/semi- 
implicit numerical scheme (ALAT; Amari et al. 1996b). We 
then start shearing the initial state, using a linear ramp of 10r A 
to reach a maximum velocity of 10” 2 at / = 10. The photo- 
spheric twisting motions are applied up to about t = 200. The 
time step used for the time advance is chosen between 0.05 
and 0.1. Thanks to our numerical scheme, this allows us to 
effect the simulation in a reasonable computational time. 

The system evolves quasi-statically through a sequence of 
force-free configurations. This type of boundary velocity profile 
allows to us build up a large amount of shear in the configu- 
ration along the neutral line (about 70°-75°, with 90° meaning 
a field aligned with neutral line), as is often observed in active 
regions (Heyvaerts & Hagyard 1 99 1 ). At / = 200, the boundary 
velocity field is progressively switched off (with a linear ramp 
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Fig. 1 . — Configuration obtained at/ = 550 after applying a localized bound- 
ary twisting vector field shown at the bottom of the figure (superposed to the 
photospheric image of the B. distribution) to the equilibrium configuration 
obtained at the end of shear buildup (shearing motions followed by a relaxation 
phase) at i = 400. 

of 10 time units) and the system is allowed to relax up to 
t ~ 400 to an accessible neighborhood equilibrium. 


3. TWIST BUILDUP BY IDEAL MHD PHOTOSPHERIC MOTIONS 

For t > 400, we now start applying a twisting boundary 
velocity field such that the free function \p in equation (1) is 
given by two Gaussian distributions located on each side of 
the neutral line of this new initial configuration, as shown in 
Figure 1 (superposed to the photospheric image of B z ) with 
r 0 = 1CT 2 . This velocity field corresponds to two parallel vor- 
tices rotating in the same direction and located symmetrically 
with respect to the origin O of the plane [their center is ( ± 1 .4, 
±0.5) with respect to O] and width equal to 0.4. The ramp 


function f which is used to smoothly switch on or off the 
velocity field, is chosen to be linear. The configuration evolves 
through a sequence of force-free equilibria with a monotoni- 
cally increasing energy. Figure l shows one of this intermediate 
equilibrium obtained at / = 550 after these twisting motions 
have been applied. We actually checked that the configuration 
relaxed toward an equilibrium by performing a viscous relax- 
ation procedure (see ALAT; Amari et al. 1996b). The config- 
uration clearly presents a twisted flux tube aligned along the 
neutral line and still confined by an overlaying arcade. The 
maximum value of |a| is 8.2 in our units. The concavity is 
directed upward in the central part of the tube, which implies 
(Amari et al. 1991) a configuration favorable to material 
support. 


4. DIFFUSIVE PHOTOSPHERIC MHD PROCESS 

Let us now consider again the sequence of equilibria obtained 
in § 3 with boundary shearing motions. Let us then apply 
v x (x y y, r) = 7\,(x, y, t) = 0 with / 0 e [0, 400] as boundary con- 
ditions for t > t 0 and apply a tangential electric field corre- 
sponding to diffusion of B z on the boundary {z = 0} only, while 
the domain {z > 0} is treated ideally (77 = 0). Let us take two 
particular cases: 

L / 0 = 400: In the first case, we choose as the initial state 
for this relaxation process the configuration obtained at the end 
of the process of shearing and ideal viscous relaxation consid- 
ered in § 2. As shown on Figure 2 (t - 480), reconnection 
takes place on the sheared field lines (mainly along the neutral 
line; left), producing a large twisted flux tube (right). Although 
the net twist is about the same order as for the configuration 
obtained in § 3 (almost 3ir) with an ideal process, the amount 
of flux from which the field lines that have the upward con- 
cavity originate is larger, leading to a much larger magnetic 
dip favorable for material support. Note that in the ideal MHD 
case considered in § 3, the twisting motions produce electric 
currents that tend to make the tube deviate from its axis, while 
we actually found that in the diffusive case the large electric 
currents (corresponding to = 36.6) are created by those 
currents originating from the highly sheared field lines and are 
strongly confined to the core of the twisted flux rope — along 
the tube axis — keeping the structure aligned with its axis. 



Fig. 2. — Evolution of the configuration (obtained at the end of shear buildup t = 400; left ) from an arcade-like topology to a twisted flux rope-like topology 
{right), after applying an electric field corresponding to diffusive photospheric boundary conditions (r = 480). 
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Fig. 3. — Configuration obtained at t = 100 by applying an electric field 
(corresponding to diffusive photospheric boundary conditions) to the config- 
uration picked up at time / = 20 of the previous shear buildup phase. The 
strong electric currents confined in the core of the tube are responsible for 
more twist and magnetic dips than the configuration obtained with the ideal 
MHD fluid motions. 

2. / 0 = 20; Let us now consider as a new initial state a 
configuration that is only slightly sheared. Figure 3 shows how 
after the same time as in case 1 a twisted flux tube aligned 
along the neutral line is produced. This flux rope is embedded 
into a confining potential arcade. The amount of twist is much 
larger than in case 1 (up to about 1 0 turns, although it is difficult 
to give a rigorous intrinsic definition of twist for an arbitrary 


three-dimensional configuration). The distribution of current 
(corresponding to \a\ m3X = 66.4) is strongly nonlinear, and 
strong electric currents flow along the tube axis and are re- 
sponsible for this higher twist and for a large number of dips 
in a flux rope having a smaller section than in case I. 

It appears clear that the diffusive MHD process discussed 
in this section is much more efficient for creating a flux rope 
having more dips (a bigger twist) than the purely ideal one. 
The configuration is confined by an overlaying arcade as is 
often observed in active regions. It is worth noticing that the 
fact that we have chosen 77 0 only on the boundary is not 

crucial. 

Note that the two classes of boundary conditions may not 
be considered equivalent since the corresponding boundary 
electric fields are of fluid-type nature (for the ideal process) 
while purely resistive for the second class of process (the elec- 
tric field cannot be written as a 1? x b fluid term, and its tan- 
gential component can be imposed a priori and depends on 
solar photospheric processes or the trace of the solar sub- 
photospheric physical processes). 

More various configurations (and boundary conditions) as 
well as their stability properties are currently under study and 
will be reported in a forthcoming paper. Note that our results 
may represent a new way for analyzing more consistent models 
in which more general MHD processes implying plasma effects 
could be taken into account in a following step. 
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Abstract. We present a method for reconstructing the coronal 
magnetic field, assumed to be in a non-linear force-free state, 
Irom its values given in the photosphere by vector magnetograph 
measurements. In this paper, that is the first of a series, we pro- 
pose a method that solves the boundary value problem set in the 
functional space of regular solutions (i.e., that do not contain 
current sheets). This is an iterative method introduced by Grad 
and Rubin. It is associated with a well-posed boundary-value 
problem. We present some results obtained with this method on 
two exact solutions of the magnetostatic equations, used as the- 
oretical magnetograms. Unlike some other extrapolations meth- 
ods, that are associated with ill-posed boundary value problems, 
our method allows extrapolation to arbitrarily large heights, w ith 
no blowing up due to the presence in these methods of an in- 
trinsic instability that makes errors growing up exponentially. 

Key words: rodynamics ( MHD) - Sun: corona - Sun: magnetic 
fields 


1. Introduction 

The magnetic field dominates most of the corona, and it is prob- 
ably the origin of a large variety of structures and phenomena, 
such as flares. Coronal Mass Ejections, prominences and coro- 
nal heating (Priest 1 982). Unfortunately the magnetic field is not 
yet observationally accessible in the tenuous and hot plasma that 
fills the corona (see Sakurai 1989, Amari & Demoulin 1992, and 
references therein). One possible familiar approach consists in 
solving the equations of a model (defined by some reasonable 
assumptions about the physical state of the corona) as a Bound- 
ary Value Problem (B VP), the boundary conditions being taken 
to be the measured values of the magnetic field in the denser and 
cooler photosphere: this is the so called Reconstruction prob- 
lem of the coronal magnetic field. Many problems have been 
encountered since the early attempts of Schmidt (1964), as the 
observational problems to get rid of the ambiguity that remains 


in the transverse component of the photospheric magnetic field 
(Amari & Demoulin 1992, McClymont et al. 1997, and refer- 
ences therein), or the problems related to the choice of boundary 
conditions that make a well set BVP (Aly 1989) 

In the simplest approximation the coronal magnetic field is 
current-free. This only requires the longitudinal component of 
the photospheric field as a boundary condition (Schmidt, 1 964), 
and the solution can becomputed using either a Green's function 
method or Laplace solver methods for the magnetic field or the 
vector potential. The mathematics of the various related BVPs 
(e.g., their welI-(or ill-) posedness properties), are also known 
(Aly 1987, Amari et al. 1998). 

In many active regions, where the magnetic configuration is 
knowm to have stored free energy, the current free assumption is 
not relevant. One can then introduced the so-called constant-a 
force-free hypothesis, which allows for the presence of elec- 
tric currents in the corona. The magnetic field is computed, for 
a given value of a, from its longitudinal component by using 
either Fourier transform (Nakagawa et al. 1973, Alissandrakis 
1981) or Green’s function (Chiu & Hilton 1977, Semel 1988) 
techniques. Other spectral methods have been recently proposed 
(Boulmezaoud et al. 1 998). It is also possible to solve it by regu- 
larizing an ill-posed BVP (in which the three components of the 
magnetic field are used) (Amari et al. 1998). However, the non- 
regularization or partial regularization of the so called Vertical 
Integration Method (VIM), leads in general to an amplification 
of errors (Wu et al. 1990 and refemcess therein, Cuperman et 
al. 1990a-b, Cuperman et al. 1991, Demoulin et al. 1992). In 
addition the total energy of the linear force-free field in an un- 
bounded domain such as the exterior of a star shaped domain is 
infinite, and is in general infinite in the case of the upper half 
space (except for some particular periodic solution satisfying 
some special conditions, Alissandrakis 1981, Aly 1992). More- 
over the electric currents are uniformly distributed, while ob- 
servations clearly show strong localized shear along the neutral 
line of many active region magnetic configurations (Hagyard 
1988, Hofmann & Kalman 1991). 
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Modeling such strong localized electric currents needs to 
assume that the coronal magnetic configuration is in a non- 
linear, force-free state. In this case one can distinguish two types 
of methods associated to different classes of BVPs, Extrapola- 
tion Methods and Reconstruction Methods. In the first class of 
methods the three components of the magnetic field are used 
as boundary conditions. The equations are thus vertically inte- 
grated step by step, from the photosphere towards the corona, 
without incorporating any type of asymptotic boundary condi- 
tions. This give rise to the VIM (Wu et al. 1990 and references 
therein, Cuperman et al. 1990a-b, 1991 , Demoulin et al. 1992). 
This method, associated to an ill-posed boundary value prob- 
lem, has not yet been proved to be convincingly regularized, still 
ending with an exponential growing of the errors with height, 
prohibiting extrapolation up to reasonable heights. The second 
class of methods considers a BVP that only requires the normal 
component of the field on the boundary (B n ) and the normal 
component of the electric current say, where B n > 0. Now the 
problem is considered in the whole domain and the solution is 
globally sought. It has been tackled by the use of iterative meth- 
ods introduced by Grad and Rubin (Grad & Rubin 1 958, Sakurai 
1981 , Sakurai et al. 1985) and by the Resistive MHD Relaxation 
Method (Mikic & McClymont 1994, Jiao et al. 1997). Roume- 
liotis (1997) presented a Relaxation Method in which the three 
components of the magnetic field are used at the photospheric 
level. Another method (see Amari & Demoulin 1992), is the 
Method of Weighted Residuals (Prid more -Brown 1981). This 
method is based on the minimization of two residuals, one asso- 
ciated with the Laplace force that has to vanish for a force-free 
magnetic field, and the other one with the difference between 
the directions of the observed transverse photospheric magnetic 
field and of the computed one. However, some aspects, such as 
the choice of test functions to be used for scalar products, as 
well as some other points concerning the definite positiveness 
of one functional to be minimized, are not yet clear. Other com- 
putational schemes such as collocation or least square methods 
have also been proposed in Amari & Demoulin ( 1 992), but they 
have not been tested so far. 

Sakurai (1981, 1985) presented a Green’s function approach 
of the Grad-Rubin formulation. Practically, the standard Green- 
Function formulation is however numerically expensive, since 
at each step of the iterative scheme one would need to com- 
pute an integral over the whole volume to get the value of the 
magnetic field at each point! An alternative approach proposed 
by Sakurai (1981, 1985) is to discretize the integral involving 
the Green’s function by introducing “finite-elemenf’-like di- 
cretization for the electric currents. The process thus consists 
in starting from an initial current-free field line, putting current 
on it, and then retracing the correct perturbed field line carrying 
the electric current just put on. In this method, the field lines 
are discretized into a finite number of nodes (which define the 
degrees of freedom of the problem) and the nodes locations then 
become the unknowns of the problem for tracing the field lines. 
The latter are determined by solving a system of nonlinear al- 
gebraic equations, whose convergence is related in some sense 


to the absolute value of o, and has not been proved to hold lor 
large values of o. 

In this paper we consider another class of Grad-Rubin Meth- 
ods that used the vector potential representation of the magnetic 
field. The paper is organized as follows. In Sect. 2 we present 
the general problem that is solved. In Sect. 3, we present the 
class of Grad-Rubin-like computational methods for solving 
the non-linear force-free case. We introduce in particular a new 
Vector Potential formulation in Sect. 4. We then present some 
results obtained with our method when applied to some par- 
ticular known solutions in Sect. 5. Sect. 6 gathers concluding 
remarks. 

It should be noted that a portion of the present study has been 
published in the proceeding of a conference (Amari et al. 1997). 

2. The problem 

The set of equations that describe the equilibrium of the coronal 
magnetic field in the half-space = {z > 0}, when plasma 
pressure and gravitational forces are neglected, are the well 
known force-free equations (Parker 1979): 

7xB = o(r)B in (1) 

V B = 0 in n, (2) 

in which a(r) as well as B are unknowns. 

The analysis of set of characteristics curves of this system, 
which is in general nonlinear (Grad and Rubin 1958, Parker 
1995), shows that this system has a mixed elliptic-hyperbolic 
type structure. This complex structure of the problem ( already 
known in fluid mechanics as the Beltrami field equations) makes 
this problem a formidable task to solve, and still makes it an open 
field of research in applied mathematics (Laurence & Avel- 
laneda 1993), even in bounded domains. Moreover the astro- 
physical constraints, as seeking a solution in a domain that may 
be unbounded as Q = {z > 0} add another non-trivial diffi- 
culty. 

This mixed nature implies the requirement of two types ot 
boundary conditions: 

- First of all, the elliptic part, resulting from the assumption 
that the RHS of Eq. ( I ) is given (the electric current), is rather 
well known, since it is nothing else than the Biot and Savart 
law, and just requires the value of B n on dQ to compute B 
in the whole domain, as expected for any elliptic problem: 

B n (3) 

where b (i is a given regular function. 

- Then from Eqs. ( 1 )-(2) one gets a hyperbolic equation for a 
(for B given): 

B.tfa(r)=0, (4) 

and therefore one may give the value of a in the part 
of dQ. where B n > 0, say: 


a\oQ+ - a (J , 


15) 
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where Qu is a given regular function. Note that this type ol 
boundary condition is sufficient if one reasonably assumes 
that every field line of the coronal magnetic field has its two 
footpoints connected to the boundary <9f>. Configurations 
having non-connected field lines (magnetic islands) would 
otherwise lead to the impossibility ot transporting informa- 
tion from the boundary <912 (Aly 1988). 

Because 12 is unbounded, one may also require the asymp- 
totic boundary condition: 

lim |B 1 - 0. (6) 

M-v* 

3. Grad-Rubin approach 

Let us follow the approach that was proposed by Grad and Ru- 
bin (1958). The previous underlying mixed elliptic-hyperbolic 
structure of the system of equations is exploited by introducing 
the following sequences of hyperbolic and elliptic linear B VPs: 


B (n, .Va (n) =0 in Q. 

(7) 

o*”^lan* = Q o . 

(8) 

and 


v X B ,n+1) = Q !n) B< n) in n, 

(9) 

V • B |n+1) = 0 in Q, 

(10) 

fll B+1, lflu = h , 

(II) 

lim |B (n+1) | = 0. 

(12) 

|rhx 


with B° the unique solution of: 


V x B n = 0 in 12, 

(13) 

(an = fc) , 

(14) 

3 

B 

ii 

o 

(15) 

1 r | ►oc 


that is given by (Aly 1989): 


W 

II 

< 

(16) 

x 1 f n , „ dx'dy' 

,,r| -iL“' lr| |r- r '| ■ 

(17) 


These sequences of problems must be proved to converge 
towards the solution of the original BVP defined by Eqs. (1) 
and (2) provided with the set of boundary conditions given by 
Eqs. (3), 5). One can use them to address theoretical issues such 
as i) existence of solution ii) uniqueness iii) continuity of the so- 
lution with the respect to the boundary conditions. These three 
points define a well-posed BVP in the sense of Hadamard ( 1 932) 
and has been discussed for other BVP associated to extrapola- 
tion methods (Low & Lou 1991, Amari et al. 1998). Note that 
the last point is important because of the presence of errors 
in the measurements of the photospheric magnetic field and of 
the possible non-force-free character of the field at the photo- 
spheric level, where pressure and dynamic forces can play a non- 
negligible role (Aly 1989, McClymont et al. 1997). Of course 
those three points depend on the functional space in which one 


seeks the solution, and in particular on the the regularity of the 
solution (Amari 1991). 

Bineau (1972), considered this BVP in the Holder func- 
tional spaces (set of functions sufficiently regular and whose 
derivatives are also regular enough, Brezis 1983). The BVP is 
then proved to be well-posed when a < a c . However, this 
proof rests on the following assumptions: (i) The domain 12 is 
bounded, (ii) The field Bu as well as B have a simple mag- 
netic topology (then they must not vanish in 12). It is however 
possible to show the existence of a solution for 12 bounded, in 
more general spaces (when (a.B) € L * x that is 

in a functional space such that solution may admit separatrices 
surfaces, null points, and current sheets, (Boulmezaoud, Amari 
& Maday, in preparation). Uniqueness of the solution has not 
yet been proved. 

4, A vector potential formulation 

4.1 . Gauge for B n fixed on <912 

To ensure that B is divergence free (Eq. (2)) we use the vector 
potential representation for B. Since in BVP (10)-(12) B u \qq 
is fixed, the vector potential A should be determined such that 

B = V x A in 12, (18) 

B„|osi = 6 d < 19 > 

This representation is not yet unique, since if A is a potential 
for B then: 

A = A + V0 (20) 

where <p is an arbitrary scalar function, is also a vector potential 
for B. Uniqueness is obtained by the choice of a particular 
gauge. There are several possible choices (Dautray & Lions 
1982), but these do notingeneral take into account Eq. (1 l)(one 
well known choice is for example V * A = 0 and A n \oQ — 0). 

Our gauge is fixed by imposing that A is the unique vector 
potential such that : 

B = V x A in 12, (21) 

V * A = 0 in 12, (22) 

V £ • A t = 0 on <912 (23) 

where the subscript t in stands for the trace (when it exists) 
of the operator or the field £ on the boundary (in particular 
in cartesian coordinates: V* ■ g = di9i + djQj on the plane 
Efc = {r • hk — constant} with (?', j, k) := (x, y, z) and ru- 
standing for the unit vector normal to the current boundary plane 
E/,.). Note also that one readily gets: 

d n A n = 0 on <912, (24) 

where V n {f) = h ■ V/. The proof that A is unique is straight- 
forward since from Eqs. (20)-(23) one gets that <p is the unique 
solution of a Laplace equation: 

A0 = 0 in 12, (25) 

4> = <po on <912 (26) 
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(note that <po is also obtained once solving a Laplace equation 
on dQ, and is also unique once do is prescribed on the border 
T of the boundary an). 

Then with this choice of Gauge, A is the unique vector 
potential that satisfies: 


B = V x A in Q, 

(27) 

<1 

> 

11 

o 

P 

(28) 

v £ • At = o on dq 

(29) 

A t = V ± \ on OQ 

(30) 

d n A fl = 0 on OQ 

(31) 

where V 1 is the operator defined on 
(i.e., V x x = x n) and where \ 

OQ such that V x ■ V = 0 
is the unique solution of 

- A s x ~ h in OQ , 

(32) 


X = 0 or d n x = 0 on T. (33) 

where A* is the Laplacian operator on an (i.e. , A,/ := V^/) 


4.2. BVP for A 

One can then rewrite BVP(10)-((12) in terms of the potential 
vector A that is then the unique solution of the following BVP 
(referred to hereafter as BVP- A): 


V x A (rl) .Vo4 n} =0 in 

ft, 

(34) 



(35) 

= c*o • 


(36) 


and 


4.3. A two- level iteration procedure 

Let us define the sequence (a t) Ji< p <p and the mononotic in- 
creasing sequence (« P )i<p<P such that: 

a.\) r (x) = u p c\o(x) for x € OQ , (44) 

■«! = r, (43) 

up — 1 , (46) 

where c is a ‘’small enough" real number, and P is a ’’large 
enough" integer. 


One can then 'generate a more general sequence of linear 


BVP for (A p l \ Qp' ) )„>i.i<p<p given by: 

V x A<, n ’.Vtt<, n) = 0 in ft, 


(47) 

(48) 

a p’ , li5!i+ = Q n,, • 


(49) 

and 

- AAj, n+11 = a (n, V x A'/” in 

n. 

(50) 

A p”' +l \ — V x \ on 3ft. 


(51) 

d, l A l , (n+v> n = 0 on dQ 


(52) 

_ lim jAg ,+I) | = 0. 


(53) 


jrj— *:>c 


One may initialize the iteration procedure for p = 0, with the 
unique solution of BVP(13)-(15) (which would be equivalent 
to choose in) - 0). A possible choice for (u p )i< P <p is for P 
given: 


- AA fn+1) = a (n) V x A (n) in Q, (37) 

A'" +1) = V x x on OQ. (38) 

3 ^l A^ 1+1, =0 on OQ (39) 

lira |A (rl+1, | = 0. (40) 

}r|— *oc 


The solution A^ n+1) of the linear elliptic mixed Dirichlet- 
Neumann BVP is in general regular (A (n_l " 15 € C 2 (Q) U 
C l (OQ) :i ). 

One can then prove that: 

Vn>l;,A fn) satisfies V - A = 0 in Q, 


(55) 

One clearly notices that for every value of p one needs to solve 
a sequence of linear BVPs for all n > 0. This corresponds to 
a progressive injection of a at the boundary which turns out to 
improve convergence of the classical Grad-Rubin scheme. 

4.4. Numerical implementation 

We have developed a code called EXTRAPOL, based on the 
method described in the previous sections. 


Proof: Applying the operator V- to both sides of Eq, (37) and 
using Eq. (2) for B n , one gets: 

- A(V • A ( ” +1) ) = V • (a (n) V x A (>1) ) (41) 

= V x A (n \Va <n) = 0 in U42) 
V • A (n+l) = 0 on <9Q, „ (43) 

Whence V ■ A( n+1) = 0 is the unique solution tending to zero 
at infinity for this BVP. Note that since the initial potential mag- 
netic field B (0) clearly satisfies V ■ J< ()) = 0 (where J stands for 
the associated electric current), this property is preserved for all 
n > 0. 


- i) The computational domain Q is supposed to be the 
bounded cubic box [0. Lx] x [0, Ly\ x [0, Lz] (instead of the 
infinite upper half space), that w p e discretize as Qf t using a 
non-uniform structured mesh for finite difference approxi- 
mation. This staggered mesh used for the the various compo- 
nents of the vector potential A, the magnetic field B and a, 
is the same as the one used in our MHD code METEOSOL 
used for three-dimensional dynamic evolution (Amari et al. 
1996). 

- ii) We use as a boundary condition for BVP (37)-(40) on the 
lateral and top boundaries of the box, B n = 0, which owing 
to our gauge choice (Eqs. (22)-23)) is equivalent to impose 
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on these boundaries: 

A t = 0 (56) 

d n A n -= 0 (57) 

This type of boundary conditions, whose aim is to mimic the 
far field behaviour at infinity (as one would expect tor the 
magnetic field in the actual infinite half-space), implies that 
the top and lateral boundaries of the box have to be chosen 
sufficiently far away from the main region of interest. This 
can be achieved at relatively low cost since our mesh is not 
uniform, and therefore large cells can be put in the far-field 
region. 

- iii) The various differential operators (Eqs. (47)-(53)) are 
then dicretized on this mesh to second order accuracy. The 
Laplacian operator (in the Dirichlet-Neuman BVP (50)- 
( 53 )) leads to a 7 diagonal sparse positive definite matrix. 
The corresponding linear system is solved by use of an it- 
erative method, in which the matrix is not stored but the 
matrix-vector product is generated explicitly by the opera- 
tor (and only one more array is stored for building a precon- 
ditioner to accelerate the convergence of the method). This 
memory space saving allows the method to be implemented 
on a workstation with reasonable central memory size, and 
not only on supercomputers. We actually run the code on 
both machine types although the results presented here cor- 
respond to runs performed on a CRAY C90 machine. 

- iv) The numerical solution of Eqs. (47)-(49) is performed by 
using a characteristics method approach, since those curves 
are the field lines. Let (X; s) be the characteristics, solution 
of 

X' = B(X), X(0) = q (58) 

for q given in (the prime symbol standing for differ- 
entiation with respect to the parameter that runs along the 
characteristics). Then for any node € Qh on which a is 
defined, one gets a h as 

= ttu(Xan+ (q/i)) (^ 9 ) 

where X# u +(q h) = X(q/i,S 0 n+) is the intersection of 
{X(q; 3 ) : s < 0} with dQ + . Since cto is known at the 
nodes that do not in general coincide with ao(Xan+(q/i)), 
an interpolation from its four nearest neighbors eventually 
gives o(q/x). We have then derived tw'o methods: a) In the 
first one, once a step is chosen for field line integration, 
one goes backwards along the characteristics using a sec- 
ond order predictor-corrector scheme. Clearly one can save 
computation time by avoiding going back up to <9Q + . This 
is achieved by marking the nodes in the domain where a 
has already been computed, and then linearly interpolating 
a from its nearest neighbors as soon as'the current node is 
surrounded by such marked nodes, b) In a second method 
(Pironneau 1 988) one avoids fixing a step by using a slightly 
less accurate scheme that consists in going backwards along 
the characteristics following the faces of each cubic cell that 
is centered on an a-node, approximating the characteristic 


curve by a polygonal line made of the segments [q fc . q^ 1 ] 
where q° - q fi and q M is the intersection of the line 
{q k - /fB(q A ')} ;/ >u with the boundary dC m of the cubic 
a-cell C m that contains q^ and q^ - //B(q fc ) (with 1 ] > 0 ). 
This method is then faster than the previous one since there is 
no step size to be fixed a priori. Unlike for the first method, 
in a non-uniform mesh, each cell is crossed in ‘one step’ 
only, which makes this method faster in the big cells re- 
gion. Despite this difference in the computational speed we 
have kept the two methods available because of their slight 
accuracy difference. 

5, Application to some know n exact force-free solutions 

We now test the scheme presented in the previous section by 
running our code EXTRAPOL on some analytical and semi- 
numerical exact solutions of the non linear force-tree equations 
Eq. (l)-(2). The boundary values of these exact solutions are 
used as simulated magnetograms. Hopefully, in these cases one 
knows the solution above in the domain too, and compare the 
reconstructed and the exact solutions (which is not the case 
for the actual corona!). There are only very few- knowm exact 
solutions of the force-free equations. Let us presents the results 
obtained with our code on tw r o cases that have been also used 
by other methods such as the VIM (Demoulin et al. 1992) for 
the first one and the Resistive Relaxation Method (Mikic & 
McClymont 1994) for the second one. Note that another class 
of related solutions that will not be tested here are those found 
by Cuperman & Ditkow'ski (1991). 


5./. The Low (1982) solution 

Our first target is the well-known solution of Low (1982) for 
which the magnetic field B is given by: 


B x 

Bo . . 
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is 

related to 


a(r) = —<j> f {r) . (63) 

We choose for the function <p 


6{r) = r 0 ct m tanh(r/r 0 ), 
which owing to Eq. (63) gives: 

/ v a m 

a * tanh(r/r 0 ) 


(64) 


( 65 ) 


We fix hereafter r () — 4 and a m =0.2 
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Our numerical box size corresponds to the choice L x — 48, 
L 2 = 48,L 2 — 48. A non-uniform mesh with 51x51x51 nodes 
was chosen as in Demoulin et ai. (1992). The analytical solution 
is then computed on the mesh and in particular the values taken 
by B z and a on the boundary provide boundary conditions for 
our force-free reconstruction procedure. We choose P — 25 in 
Eq. (53) (i.e., the parameter necessary to fix the outer iteration 
corresponding to the injection of ao). Our method converges up 
to a Lorentz force of order 10 for a number of inner Grad- 
Rubin iterations N iterations — 6. The numerical error is defined 
as in Amari et al. (1998). We also found that choosing P — 15 
implies increasing Ngradrub up to about 12 to reach a Lorentz 
force of the same order. Fig. 1 shows some field lines of the exact 
solution (top) and the the corresponding field lines obtained 
from our computation. 

Some discrepancies (up to few percents) between the ex- 
act and the computed solution are found in the domain, and 
these can reach almost .2 for the field lines approaching the 
lateral boundaries of the box. These can be explained by our 
choice for the boundary condition (B n — 0) on these bound- 
aries for the computed solution, while the exact solution does 
not decrease fast enough and even more pathologically in the 
horizontal plane (see Amari et al. 1998). Note that because ap- 
plying this boundary condition results in a difference between 
the computed and exact solution, but still allows to reach a force- 
free equilibrium. However this equilibrium shows a different 
behaviour than the exact solution near the boundary, but there is 
no intrinsic instability as in the VIM (Wu 1990 and references 
therein, Cuperman et al. 1990a-b, Demoulin et al. 1 992). It is 
worth noting that we have also performed some higher resolu- 
tion run, with Nx = 101 , Ny = 101 , Nz - 101 which, unlike 
the VIM, gave even better results, allowing the boundaries to be 
pushed far away. Note that this ''robustness *' property (good 
behaviour while increasing spatial resolution) as well as the 
convergence of the method even for this type of lateral and top 
boundary conditions results from the well-posed formulation we 
have adopted, unlike for the VIM which is associated to an ill- 
posed mathematical problem (see Low Sc Lou 1991 and Amari 
et al. 1998). In this latter method errors increase exponentially 
with height (Demoulin et al. 1992) and this is a property intrin- 
sic to the method (and not the numerical scheme used for the 
extrapolation), which implies that the computed solution will 
eventually diverge, while our solution never diverges for an ar- 
bitrarily large box. Actually the bigger is our box, the bigger the 
region of agreement between our solution and the exact one is, a 
property that we checked with the higher resolution run, pushing 
the lateral boundaries to L x — 120, L z = 120, L z = 120. 

5.2. Low & Lou 's (1991 ) solution 

We have also tried the particular case of the exact force-free 
solution presented in Low & Lou (1991). Unlike Low's (1982) 
solution, it requires some numerical calculations. 

The solution is supposed to be axi-symmetric and writes in 
spherical coordinates: 



Fig. 1. Example of force-free reconstruction with the Vector Potential 
Grad Rubin Method compared to the exact analytical Low’s (1982) 
solution. The computed solution {bottom), matches the exact solution 
{top) up to few percents in most of the interior of the computational 
box. Some discrepancies occur for field lines near the lateral and top 
boundaries because of the slow asymptotic decreasing behaviour of 
this particular solution, while the boundary condition B n = 0 has 
been imposed on the boundaries in the computation. 
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( 68 ) 


where Q{A) is an a priori unknown function of A(r. Q), a so- 
lution of the nonlinear partial differential equation (see Low & 
Lou 1991). A family of solutions can be generated by choosing 


P(cos0) 


— n a *+ « 


Q — aA 


(69) 

(70) 


for odd n, and a a real constant. P is then the solution of the 
following boundary -value problem: 

d 2 P 

{l ~ cos 29 ) d(^ +n{n + l)P 

= o. (71) 

n 

P(-l) = P(l) = 0. (72) 



We then solve numerically Eqs. (71)- (72). Usual transforma- 
tions (Low & Lou 1991) then allow to get the solution in carte- 
sian coordinates, in the upper half space. 

Our numerical box is taken such that L x = 8, L z — 8, 
L z = 4. A non-uniform mesh is generated with N x = 60, 
N y — 60, N z = 40 with most of the cells concentrated in the 
inner stronger field region. Once BVP (7l)-(72) is solved, one 
deduces the corresponding three components (B Xy B y ,B z \ the 
associated electric currents and a = ^ on the same nodes 
(xh, yh , Zh) of the mesh used by our force-free reconstruction 
code EXTR APOL, and then computes the solution. One then use 
B z {xh, yh , 0) and a(x h, 0) (for the nodes in dQ + only) as 

boundary condition for the reconstruction procedure. We found 
that using P = 20 and 4 inner iterations( AT iterQtiorij5 — 4) 
allows to decrease the Lorentz force down to values of order 
lCrT 

Fig. 2 shows some field lines of the exact solution (top) and 
the corresponding field lines resulting from our reconstruction 
procedure. The errors, defined as for the previous case (Low’s 
( 1 982) solution), are even less or of the order of 1% in the larger 
part of the domain, except again near the lateral and top bound- 
aries where the imposed boundary condition B n = 0 and the 
exact one disagree. Actually, those discrepancies are however 
smaller than those of the case of Low's ( 1 982) model for the lat- 
eral boundaries because the magnetic field now decreases faster 
with distance. The case of the top boundary is different because 
of the existence, in the exact solution, of a pathological field- 
line in the center of the box that crosses almost vertically the 
top boundary while it has to match the applied boundary condi- 
tion B n - 0 in the calculation, which will be difficult to fulfill, 
even with a large box. Note that despite the much better asymp- 
totic behaviour of this force-free solution for the magnetic field 
the electric currents are distributed on a scale that is still large, 
which results in a configuration that does not quickly approach 



Fig, 2. Non linear force-free reconstruction ( with the Vector Potential 
Grad Rubin Method) of the semi-numerical exact solution of Low & 
Lou (1990). The computed solution ( bottom ), and the exact solution 
(top) agree in most of the computed area. The existence of a pathologi- 
cal field line (in the exact solution) that crosses almost perpendicularly 
the top boundary, implies larger errors near this boundary since the 
computed solution corresponds to the boundary condition B n = 0. 
The boundaries of the box are put far away enough from the inner 
stronger field area. 
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toward the potential field as it is often the case in the corona, 
outside regions of more localized electric currents. 

6. Conclusions 

In this paper, we have presented a numerical method for recon- 
structing the coronal magnetic field as a force-free magnetic 
field from its value given on the boundary of the domain. Let 
us summarize here the main points we have discussed and our 
main results; 

(a) The boundary-value problem is formulated such that it cor- 
responds to a well-posed mathematical problem: the normal 
component of the magnetic field is imposed on the boundary 
of the domain, and a only on that part of this boundary where 
B n > 0. We impose B n = 0 on the lateral boundaries so that a 
does not need to be specified on these boundaries, provided that 
these boundaries are put far enough to mimic the behaviour of 
the solution at infinity. 

(b) We have derived a Grad Rubin Vector Potential formulation 
of this B VP to ensure div B = 0 up to machine roundoff numer- 
ical errors. We have shown that this problem may be equivalent 
to solve a sequence of linear elliptic boundary-value problems 
for the Vector Potential, and hyperbolic ones for a. This cur- 
rent formulation is relevant for seeking regular enough solutions 
but not equilibria having current sheets. Weak formulations of 
these methods are however currently under study and would 
be reported in a next Paper of the series (Amari & Boulmeza- 
oud 1999, in preparation). We have implemented this method 
in our computational code EXTRAPOL, in a relatively efficient 
numerical way. Other mathematical approaches allowing the ex- 
istence of critical points in the configuration are also currently 
studied. 

(c) We have successively applied our method to theoretical mag- 
netograms obtained from two exact known solutions, the so- 
lutions of Low (1982) and Low & Lou (1991). The method 
converges up to a small residual Lorentz force, in a reasonable 
number of iterations. Some discrepancies between the exact so- 
lution and the reconstructed one occurred near the top or lateral 
boundaries of the computational box, and have been explained 
by the relatively bad asymptotic behaviour of Low's (1982) so- 
lution, or the existence of an almost vertical pathological field 
line in the solution of Low & Lou (1991), which makes difficult 
to match our applied boundary conditions on these boundaries 
( B n — 0). Other approaches involving the assumption of po- 
tential field near the boundaries, or approximation of the Green 
formulae that can explicitly give the normal component at those 
boundaries are currently underdevelopment. Another approach 
could be to map the infinite upper half-plane onto the bounded 
square box by using a class of mappings that represent the gen- 
eralization of conformal mappings used irrtwo dimensions. 

(d) Our formulation is better than the (VIM) (Wu 1985, Cuper- 
man et al. I990a-b, Demoulin et al. 1 992) since it corresponds 
to a mathematically well-posed boundary value problem. Al- 
though it may exhibit some residual discrepancies with the ex- 
act known solution, errors never increase exponentially up to 


blowing up as in the VIM. Moreover as it was shown by Bineau 
(1972) another consequence is that the solution is expected to 
be continuous respect to boundary conditions, at least for a not 
to large (Amari et al. 1998). 

(e) Our method is different from Sakurai’s (1981) approach in 
which, instead of solving an elliptic problem for B, he uses 
a more local approach where the location of the nodes that 
discretize a given field line are computed once some electric 
currents (a) are injected in this field line, as the solution of 
non-linear system of equations that does not take into account 
the contribution, of the whole computational domain (as one 
would expect in an elliptic problem). This approach allows a fast 
enough computation, which might be useful for some very con- 
centrated (almost thin isolated) tube-like configurations, but it is 
not yet clear how this truncation procedure (by solving a single 
problem for each field line) may be involved in the numerical in- 
stabilities encountered in solving the nonlinear system forcases 
corresponding to large values of a. Indeed Sakurai’s (1981) 
approach might be considered as a Lagrangian discretization 
method while we have presented a Eulerian type discretization 
that would be more suited to highly sheared magnetic configura- 
tions. The two methods should be worth to be kept and used for 
different types of data and configurations. The results presented 
in this Paper seem to be optimistic as regards the application of 
the method to simulated magnetograms. The next step currently 
underdevelopment is the application of this method to various 
sets of data provided by vector magnetographs. However there 
are several important points that need to be emphasized, and 
that make actual data much more difficult to handle than exact 
force-free solutions: 

i) First of all data are much more noisy, because of the errors on 
the transverse magnetic field measurements that are larger than 
on the longitudinal one (Amari & Demoulin 1 992, Klimchuk & 
Canfield 1993, McCIymont et al. 1997). Other errors may also 
arise after the resolution of the 180° ambiguity that exists on 
the transverse component. These errors depend on the method 
that is used (Mikic & Amari 1999, in preparation). Eventually 
the non- force- free character of the photosphere ( Aly 1 989) may 
be taken into account. Actually from point (b) above, the well- 
posedness of our formulation (for at least a not too large), would 
make the solution not very sensitive to errors expected on the 
photospheric measurements. We are currently working on the 
project of simulating the error effects (Amari et al. 1999, in 
preparation) of these instrumental errors on the transverse field 
components, by introducing some random noise in the simulated 
data obtained from some highly sheared force-free solutions 
obtained by a relaxation code (Yang et al. 1986, Klimchuk & 
Sturrock 1992), and then reconstructing them w-ith our method. 

ii) Related to these errors, one may also find that, unlike theoreti- 
cal magnetograms, actual data are far from smooth. This implies 
that any reconstruction method should be either robust or one 
will have to smooth the data prior to reconstruction, which may 
introduce possible added deviations from the sought solution, 
since there is no unique way of smoothing. 

iii) One non negligible difficulty that has to be taken into account 
is the needs for computing a from the photospheric normal 
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components of the magnetic field and of the electric current. 
Weak field regions cannot be ruled out in a straightforward way 
since high shear can be localized near the neutral line (Hagyard 
1988). " 

iv) One final point is that unlike for theoretical magnetograms, 
one never knows a priori the solution in the corona in order 
to check the reconstructed one. However an alternative can be 
the use ofYOHKOH or SOHO/EIT data (for different heights). 
These data would have to be used a posteriori to check if the 
computed structures has such loops or ‘’footpoints” that match 
the coronal observed ones, but not use these data set to fix a 
remaining free parameter such as a in linear force-free constant- 
a extrapolations! 
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